Determination of Characteristics and Establishment of Discrete Element Model for Whole Rice Plant

Abstract

In order to accurately establish a discrete element model for the whole plant flexibility of upright rice during the harvesting period, several physical characteristics, such as geometric features, moisture content, and density, of the entire rice plant were measured, along with frictional properties, such as the static and rolling friction coefficients, and mechanical properties, including the elastic modulus and restitution coefficient. A flexible and upright discrete element model of the rice plant was established using the DEM method based on the Hertz–Mindlin (no slip) and Hertz–Mindlin with bonding mechanical models. The parameters were optimized through Plackett–Burman screening experiments, steepest ascent experiments, and Box–Behnken optimization experiments to accurately determine the discrete element model parameters of each component of the rice plant. The calibration process of the contact parameters between rice grains and steel was analyzed in detail as an example, resulting in a calibration error of 0.68% for the natural repose angle. Taking the calibration of the contact parameters between the main stem and steel as an example, a detailed analysis of the calibration process was conducted. The calibration resulted in a calibration error of 2.76% for the natural repose angle and 2.33% for deflection. This study lays the foundation for understanding the mechanical response of rice and machinery when they are coupled together. Additionally, it provides valuable references for establishing discrete element models of plant species other than rice.

1. Introduction

Rice is one of the four major staple crops in China, with a planting area of 30 million hectares throughout the year. It is also one of the most important food crops in the world [1,2]. Rice production involves various agricultural processes, such as planting, management, harvesting, and storage [3,4]. During the harvesting stage, the interaction between combine harvesters and rice plants is necessary to obtain clean rice grains [5,6]. Because of the flexible and upright characteristics of whole rice plants during harvesting, studying the characteristics of whole rice plants provides important guidance for the development and design of related harvesting machinery.

The discrete element method (DEM) is widely used in the field of agricultural engineering to explore the interaction between plants and machinery [7,8,9]. The accurate establishment of a plant’s discrete element model is a prerequisite for in-depth analysis and exploration of related mechanisms. Sadrmanesh et al. [10] employed the DEM to develop a numerical model to simulate the tensile behavior of plant fibers, where the fiber model was constructed by combining spherical particles. Guo et al. [11] combined discrete element analysis with physical experimental tests to establish a discrete element mechanical model of banana stems, providing new insights into the microscopic mechanics and nonlinear damage behavior of banana stems. Research on rice has mainly focused on a specific component, with most studies focusing on the characteristics of rice grains [12,13,14]. Liu et al. [15] accurately established a discrete element model of rice grains using three-dimensional laser scanning and particle filling. Wang et al. [16] explored the dynamic characteristics and steady-state features of rice grains with different moisture contents under stacked conditions using the discrete element method. As rice grains have smooth surfaces and lack flexibility [17], the aforementioned studies adopted the Hertz–Mindlin (no slip) model during the modeling process. Moreover, when studying the interaction between rice plants and harvesters during the harvesting stage, the focus was mostly on the threshing and cleaning stages. In studies focusing on the threshing stage, it was common to model rice grains and stems separately to investigate their behavior inside the threshing drum [18,19,20,21]. In the cleaning stage, the discrete element method was mainly used to simulate the quantitative dispersion and migration characteristics of rice grains and short stems on a vibrating screen [22,23,24,25]. The straw models in these studies were built using the Hertz–Mindlin (no slip) model, which does not account for the flexibility of different parts of rice plants. Some scholars have attempted to incorporate the flexible characteristics of straw in various ways and achieved certain results. Liu et al. [26] established a flexible straw model based on the Hertz–Mindlin with bonding model, where the proposed model was constructed by connecting straw elements (filled using the multiball method) through parallel bonds. Lenaerts et al. [27] employed a straw discrete element model composed of connected spheres, which allowed for bending. However, these studies only dealt with straws with the same physical properties, and few researchers have achieved complete modeling considering the significant physical property differences among different components of whole rice plants.

In summary, the discrete element method is commonly used to establish individual models for rice grains or straw, while fewer studies have focused on establishing a discrete element model for the entire rice plant. This limitation severely hampers the visualization and simulation of the interaction between rice and machinery during the harvesting process. The main difficulties in establishing a discrete element model for the entire rice plant are (1) how to connect various components of the rice plant according to their corresponding coordinates to form a cohesive whole that resembles the real rice plant and (2) how to arrange multiple rice plants in an orderly and upright manner on the ground, consistent with practical operating conditions.

In this study, precise measurements of the physical properties, friction characteristics, and mechanical properties of the entire rice plant were carried out. Based on these measurements, a three-dimensional coordinate system and application programming interface (API) method were utilized to establish a discrete element model for the flexible upright rice plant. Calibration of the grain, main stem, and branch parameters of the entire rice plant was conducted separately using bench-scale experiments and the discrete element method.

2. Materials and Methods

2.1. Determination of Whole Plant Rice Parameters during the Harvest Period

2.1.1. Determination of Physical Properties

(1)

Typical Geometric Characteristics of Whole Rice Plants

The “Longjing 29” rice variety, widely cultivated in Northeast China, was chosen as the experimental sample. This variety possesses the basic characteristics of conventional rice and is representative. The entire rice plant is mainly composed of grains, main stems, stem leaves, primary stems, and secondary stems, as shown in Figure 1a. Leaves do not affect the kinematic and dynamic properties of rice during the harvesting process; therefore, their presence was ignored in the measurement and modeling processes in this study.

A vernier caliper (Ningbo Deli Co., Ltd., Ningbo, China) was used to measure the dimensions of the main stems, primary stems, secondary stems, and grains. Ten samples were taken for each component, and SPSS software (version 22.0, IBM Corp., Chicago, IL, USA) was used to calculate the overall mean and standard deviation. The influence of moisture content on the size distribution of each component was disregarded during the measurement process.

Taking the centroid of the rice grains as the origin (O) of the coordinate system, the width direction (W) of the rice grains was defined as the x-axis, the thickness direction (T) was defined as the y-axis, and the length direction (L) was defined as the z-axis. Thus, a three-axis spatial coordinate system O-XYZ was established to measure the overall dimensions of the rice grains, as shown in Figure 1b. As the main stem has a hollow cylindrical shape, the length, inner diameter, and outer diameter of the main stem were measured separately. The length and outer diameter of the primary and secondary stems were also measured separately, as shown in Figure 1c,d.

(2)

Moisture Content

The moisture content of the whole rice plant directly affects its frictional and mechanical properties [28]. Because of significant variations in the moisture content among different parts of rice during the harvesting period, it is necessary to measure the moisture content of each component of rice separately. The results are presented in the form of the minimum, maximum, average, and standard deviation values. Additionally, apart from geometric and functional characteristics, the physical properties of primary and secondary stems are identical. Therefore, in the subsequent research, the primary and secondary stems were merged into a single category referred to as branches.

The drying method [29] was employed to determine the moisture content of each component of the whole rice plant. The test materials included an electric hot-air drying oven (Model: DZF−6051, Qingdao Mingbo Environmental Protection Co., Ltd., Qingdao, China) and an electronic analytical balance (Model: FA1004, Shanghai Hu Jiao Stationery Co., Ltd., Shanghai, China).

(3)

Density

The density of each component of rice during the harvesting period can be measured using the immersion method [30]. The results are presented in the form of the minimum, maximum, average, and standard deviation values. During measurement, as the densities of the main stems and branches are lower than that of water, it is required to use an object with a density higher than that of water to help submerge them in water. Eventually, the volume of the auxiliary object should be subtracted.

2.1.2. Determination of Friction Characteristics

Based on the contact materials with rice during the harvesting process, this study selected a steel plate, rice grain board (where rice grains were neatly arranged and adhered using tweezers to minimize gaps), main stem board, and branch board for the test, as shown in Figure 2a–d.

(1)

Static Friction Coefficient

The static friction coefficient was determined as shown in Figure 2e. The results are presented in terms of the minimum value, maximum value, average value, and standard deviation.

(2)

Rolling Friction Coefficient

The determination of the rolling friction coefficient is represented in terms of the minimum value, maximum value, average value, and standard deviation. The experimental materials include a self-built friction coefficient testing platform and a high-speed PCO camera (PCO Company, Bavaria, Germany), as shown in Figure 2f. Based on the principle of mirror reflection imaging, a spatial reference wall and a reflecting mirror surface were designed, and a three-dimensional coordinate system O-XYZ was simulated and established. The high-speed camera was placed directly in front of the spatial reference wall, and by measuring the spatial displacement of the material within the spatial reference wall and the reflecting mirror surface, the rolling friction coefficient of each component of the rice could be analyzed and calculated.

In the state of pure rolling at the critical point, the gravitational force along the slope direction combines with the static friction force to generate positive work on the rice grains. As the rolling action progresses, this work gradually transforms into translational kinetic energy and rolling kinetic energy of the rice grains. At the same time, the rolling friction torque performs negative work, converting the gravitational potential energy into translational kinetic energy, rolling kinetic energy, and work performed to overcome the friction, as follows:

$$\{\begin{array}{l}U=E_k+E_g+W_q\\ E_k=\frac{1}{2}m{v}_t{}^2\\ E_g=\frac{1}{2}I{\omega }^2\\ W_q=\mu mgL\cos \theta \end{array}$$

(1)

where U is the gravitational potential energy of the rice grains, J; E_k is the translational kinetic energy of the rice grains, J; E_g is the rolling kinetic energy of the rice grains, J; W_q is the work performed to overcome the rolling friction torque, J; m is the mass of the rice grains, g; v_t is the instantaneous velocity of the rice grain center point, m/s; I is the moment of the inertia of the rice grains, g∙m²; ω is the rotation angular velocity of the rice grains, rad/s; L is the rolling displacement of the rice grains along an inclined plane, m; θ is the angle of inclination, (°); g is the gravitational acceleration, m/s²; μ is the coefficient of the rolling friction between the rice grains and contact materials.

If the measured material is the main stem of the rice, it can be equivalently treated as a hollow cylinder. Therefore, the rotational inertia of the main stem of rice is:

$$I=\frac{1}{2}m(R_1^2+R_2^2)$$

(2)

where R₁ is the distance from the rotation center of the main stem to the friction surface, m; R₂ is the distance from the rotation center of the main stem to the inner epidermis, m.

Combing Equations (1) and (2) obtains:

$$\mu =\tan \theta -\frac{v_t^2(2R_2^2+R_1^2)}{2gL{R}_2^2\cos \theta }$$

(3)

Taking the analysis of the measured coefficient of the rolling friction of rice grains as an example, deviations from linear motion may occur. Selecting a stable rolling state of the rice grain and recording its initial position coordinates, the coordinates on the spatial reference wall are (x₀, z₀), and the coordinates on the reflective surface are (y₀, z₀). When the rice grain rolls steadily for n frames, the coordinates of its center point on the spatial reference wall are (x_n, z_n), and on the reflective surface are (y_n, z_n), converting the actual distance between two points in space through proportion as:

$$L=k\cdot \sqrt{{(x_n-x_0)}^2+{(y_n-y_0)}^2+{(z_n-z_0)}^2}$$

(4)

where k is the ratio of the actual size to image size.

To investigate the instantaneous velocity (v_t) of the rice grain’s center point in the n-th frame, the coordinates of the center point in the (n − 1)-th frame on the spatial reference wall are recorded as (y_n−1, z_n−1) and on the reflective surface as (x_n−1, z_n−1). Also, the coordinates of the center point are recorded in the (n + 1)-th frame on the spatial reference wall as (y_n+1, z_n+1) and on the reflective surface as (x_n+1, z_n+1). Given that the time interval between the two frames is $\Delta t$, the absolute value of the instantaneous velocity v_t of the rice grains center-point can be expressed as follows:

$$v_t=\frac{k\cdot (\sqrt{{(x_n-x_{n-1})}^2+{(y_n-y_{n-1})}^2+{(z_n-z_{n-1})}^2}+\sqrt{{(x_{n+1}-x_n)}^2+{(y_{n+1}-y_n)}^2+{(z_{n+1}-z_n)}^2})}{2\Delta t}$$

(5)

By collecting the coordinates of the center points between different frame images, we can substitute Equations (4) and (5) into Equation (3) to calculate the rolling friction coefficient between rice grains and various materials.

2.1.3. Determination of Mechanical Properties

By combining various experimental methods, such as compression tests, collision tests, and shear tests, and utilizing high-speed photography technology, measurements of the elastic modulus and restitution coefficient of each component of rice were conducted. With the help of a self-built experimental device, the mechanical damage limit conditions and mechanical properties of each component of the whole rice plant were explored.

(1)

Elastic Modulus

Because of the physical differences among rice grains, main stems, and branches, there are certain variations in the methods used to determine and calculate the elastic modulus.

The experimental instrument used for this purpose was the TA XT Plus C Texture Analyzer (Micro System Company, London, UK), as shown in Figure 3a. Rice grains possess long and short axes, and to ensure consistency during the test and calibration processes, this study measured the elastic modulus of rice grains in their natural state (with the short axis perpendicular to the horizontal plane).

According to the requirements of the ASAE S3684 DEC2000 (R2008) standard, the Hertz formula was used to calculate the elastic modulus E of rice grains:

$$E=\frac{0.338K_V^{3/2}F(1-v^2)}{D^{3/2}}{(\frac{1}{R_V}-\frac{1}{R_V^{\prime }})}^{1/2}$$

(6)

where F is the load, N; D is the deformation of rice grains, mm; v is the Poisson’s ratio of rice grains; R_V and $R_V^{\prime }$ are the curvature radius at the contact point on the upper surface of rice grains during compression, mm; rice grains in the thickness direction are subjected to squeezing by two rigid plates, and the main radius of curvature at the contact point can be expressed as: $R_V=(W^2/4+T^2)/2T$ and $R_V^{\prime }=(L^2/4+T^2)/2T$; L, W, and T represent the length, width, and thickness, respectively, of the three-axis dimensions of rice grains, mm; K_V is the coefficient determined by the radius of the principal curvature.

According to Hertz’s contact theory:

$$\cos \theta =\frac{1/R_V-1/R_V^{\prime }}{1/R_V+1/R_V^{\prime }}$$

(7)

By substituting the main radius of the curvature, we can obtain the following equation:

$$\cos \theta =\frac{L^2-W^2}{L^2+W^2+8T^2}$$

(8)

where $\theta$ is the angle between the main plane at the contact point of the upper surface of the rice grain and the rigid plate head (°).

Based on the Hertz contact theory, when the convex surface of a rice grain comes into contact with an infinitely large plane, it forms an elliptical contact area. The compressive stress is determined by calculating the maximum contact stress at the center of the contact area, which is 1.5 times the average contact stress. This calculation can be carried out using Equation (9).

$${\sigma }_{\max }=\frac{1.5F}{\pi ab}$$

(9)

where ${\sigma }_{\max }$ is the compressive stress of rice grains, MPa; a and b are the length of the major and minor axes of the contact ellipse, respectively, mm.

Equation (10) can be used to determine the contact between rice grains and the steel pressing head.

$$\{\begin{array}{l}a=m{[\frac{3F(1-{\mu }^2)}{2E}{(\frac{1}{R_V}+\frac{1}{R_V^{\prime }})}^{-1}]}^{1/3}\\ b=n{[\frac{3F(1-{\mu }^2)}{2E}{(\frac{1}{R_V}+\frac{1}{R_V^{\prime }})}^{-1}]}^{1/3}\end{array}$$

(10)

where m and n are coefficients related to the eccentricity of an ellipse, they can be obtained in the ASAE S368.4 DEC2000 (R2008) standard. The data not listed in the table are calculated using interpolation.

The main stems and branches of rice plants exhibit flexible characteristics, primarily bending during the harvesting process. The elastic modulus of the main stems and branches of rice plants was determined using a three-point bending method.

(2)

Restitution Coefficient

The restitution coefficient was determined and represented in the form of the minimum value, maximum value, mean value, and standard deviation. According to the principle of mirror reflection imaging, a spatial reference wall and a reflective mirror were positioned at a 135° angle. A coordinate grid paper was pasted on the spatial reference wall to facilitate the calibration and measurement of the data captured by high-speed cameras during the collision and bouncing of rice grains, as shown in Figure 3b.

Assuming that the rice grain impacts and falls from the adhesive plate with a certain initial velocity, it collides with a horizontal material surface. A horizontal spatial cartesian coordinate o-xyz was established, as shown in Figure 3c. Before the collision, the rice grain undergoes uniform accelerated motion in the vertical direction.

The coordinates of the rice grain in the spatial reference wall before the collision, denoted as (y_m, z_m), correspond to the coordinates in the reflective mirror, represented as (x_m, z_m). At the moment of collision, the coordinates of the rice grain in the spatial reference wall are (y₀, z₀), which correspond to the coordinates in the reflective mirror as (x₀, z₀). At this point, x₀ is equal to x_m, and y₀ is equal to y_m.

The coefficient of restitution between rice grains and horizontal materials is:

$$e=\frac{v_t}{v_0}=\frac{2t_m\sqrt{v_x^2+v_y^2+v_z^2}}{2(z_m-z_0)+g{t}_m^2}$$

(11)

where v_x represents the post-collision velocity of the rice grain in the x-direction, m/s; v_y represents the post-collision velocity of the rice grain in the y-direction, m/s; v_z represents the post-collision velocity of the rice grain in the z-direction, m/s; $t_m$ represents the time elapsed before the collision, s; $v_t$ represents the separation velocity of the rice grain after the collision, m/s; $v_0$ represents the instantaneous velocity of the rice grain before the collision, m/s.

Because of the higher moisture content and larger mass of the main stems and branches during the rice harvesting process, the coefficient of restitution between the two materials colliding on the horizontal material surface was determined by the free-fall method. The experimental device still employed the test platform for measuring the coefficient of restitution of rice grains, with the impact device removed.

The restitution coefficient between the main stems or branches and the horizontal material is:

$$e=\frac{v_t}{v_0}=\frac{\sqrt{v_x^2+v_y^2+v_z^2}}{\sqrt{2gH}}$$

(12)

The measurement methods, which are not elaborated in detail in this manuscript, can be found in

Table A1. Measurement method for each component of the whole rice plant.

Measurement Items	Determination Method	Calculation
Moisture content	(1) Place rice grains (main stem or branch) into an aluminum box and record the initial mass as Mai (i = 1, 2, 3).   (2) Put the box into an electric hot air-drying oven and dry at a temperature of 105 °C for 3–4 h.   (3) Take out the aluminum box and cool it to room temperature, then record the final dried mass as Mbi (i = 1, 2, 3).	$C_i=\frac{M_{ai}-M_{bi}}{M_{ai}}\times 100\%$ where Ci is the moisture content of each component, %; Mai is the mass of materials before drying, g; Mbi is mass of materials after drying, g; i = 1 is rice grains, i = 2 is main stem, and i = 3 is branches.
Density	(1) Weigh a certain amount of rice grains (main stem or branch) using an electronic analytical balance, and record the mass as Mi (i = 1, 2, 3);   (2) Measure a volume of pure water, Va, using a volumetric cylinder.   (3) Gently place the same amount of rice grains (main stem or branch) used in step 1 into the volumetric cylinder and record the reading of the water level as Vbi (i = 1, 2, 3).   (4) Calculate the actual volume of the rice grains (main stem or branch) by taking the difference between the volume of pure water in the volumetric cylinder, Vbi (i = 1, 2, 3) and the volume of pure water used in step 2, Va.   Note: during measurement, since the density of the main stems or branches is smaller than that of the water, it is necessary to use an item with a density greater than that of water to assist in immersion. The final volume should subtract the volume of the assisting item.	${\rho }_i=\frac{M_i}{V_{bi}-V_a}$ where ${\rho }_i$ density of each component, g/cm3; $M_i$ is the mass of a certain amount of material, g; $V_a$ is the volume of purified water, cm3; $V_{bi}$ is the volume of purified water after adding materials, cm; i = 1 are the rice grains, i = 2 are the main stems, i = 3 are the branches.
Static friction coefficient	(1) Use double-sided tape to adhere the test materials (rice grains, main stem, and branch) onto the surface of material boards, creating a rice grain board, a main stem board, and a branch board. Create steel plates with the same material as the key components that come into contact with rice during harvesting and replace the material boards after each set of experiments.   (2) Adjust the height adjustment column to ensure that the material boards are in a horizontal position, and place the test materials on the material boards.   (3) Adjust the height adjustment column at a speed of 2 mm/s to create a relative sliding tendency for the test materials. Record the angle θ between the material boards and the horizontal plane at this point.	$\mu =\tan \theta$ where $\mu$ is the static friction coefficient between various components of the entire rice plant and the contacting material; $\theta$ is the angle between the material surface and the horizontal plane when the rice plant exhibits a tendency of relative sliding on the material surface, (°).
Rolling friction coefficient	(1) Place the material board and adjust the height adjustment column to tilt the material board at a certain angle until the test material can roll on the material board. Record the angle between the material board and the horizontal plane at this point.   (2) Turn on the high-speed camera and set the frame rate to 800 frames/s. Capture images with a resolution of 1296 × 1024 pixels, with pixel dimensions of 11 μm × 11 μm. Set the exposure time to 1.25 ms, ISO sensitivity to automatic, and shutter mode to global shutter.   (3) Place the test material stationary on the material board and simultaneously start the high-speed camera to capture the rolling test material. Store the captured images in .jpg format on the computer.   (4) By calculating the ratio of the center displacement difference to time, the change in the velocity can be determined.	$\mu =\tan \theta -\frac{v_t^2(2R_2^2+R_1^2)}{2gL{R}_2^2\cos \theta }$ where R1 is distance from the rotation center of the main stem to the friction surface, m; R2 is the distance from the rotation center of the main stem to the inner epidermis, m; L is the rolling displacement of rice grains along an inclined plane, m; θ is the angle of inclination, (°); g is the gravitational acceleration, m/s2; μ is the coefficient of rolling friction between rice grains and contact materials.
Elastic modulus (rice grains)	(1) Set the parameters of the texture analyzer experiment: adjust the texture analyzer to compression mode, use a flat compressive head, and set the compression speed to 0.3 mm/s and the maximum loading distance to 1.5 mm.   (2) Place the rice grains naturally on the base, start the texture analyzer, and under the action of the compressive head, the rice grains will bear pressure and undergo deformation.   (3) The computer automatically plots the curve of load vs. displacement. Since the size of rice grains is not exactly the same, when the curve undergoes significant sudden changes, the loading should be manually stopped.	$E=\frac{0.338K_V^{3/2}F(1-v^2)}{D^{3/2}}{(\frac{1}{R_V}-\frac{1}{R_V^{\prime }})}^{1/2}$ where F is the load, N; D is the deformation of the rice grains, mm; v is the Poisson’s ratio of the rice grains; RV and $R_V^{\prime }$ are the curvature radius at the contact point on the upper surface of the rice grains during compression, mm; rice grains in the thickness direction are subjected to squeezing by two rigid plates, and the main radius of curvature at the contact point can be expressed as: $R_V=(W^2/4+T^2)/2T$ and $R_V^{\prime }=(L^2/4+T^2)/2T$; L, W, and T represent the length, width, and thickness, respectively, of the three-axis dimensions of rice grains, mm; KV is the coefficient determined by the radius of the principal curvature.
Elastic modulus (main stem and branches)	(1) Set the parameters of the texture analyzer experiment: adjust the texture analyzer to compression mode, use a flat compressive head, and set the compression speed to 0.3 mm/s and the maximum loading distance to 3 mm.   (2) Prepare the rice main stem and branches: cut the main stem or branches into sections with a length of 40 mm, and place them on the bottom support.   (3) Start the texture analyzer, and under the action of the compressive head, the main stem or branches will bear pressure and undergo deformation.   (4) The computer automatically plots the curve of load vs. displacement. When the curve undergoes significant sudden changes, the loading should be manually stopped.	$E_a=\frac{F{S}^3}{48WI}$ where Ea is the elastic modulus of the stem bending, MPa; F is the loading force, N; S is the length between two supports, mm; W is the bending deflection of the stem, mm; I is the moment of inertia of the stem cross-section relative to the neutral axis, mm4.

Following the above steps, repeat with the sample for each component ten times to calculate the measurement indicators. Except for the size and moisture content, the results are expressed in the form of the minimum value, maximum value, average value, and standard deviation.

. The corresponding measurement results are presented in

Table A2. Geometric characteristics and dimensions of the whole rice plant and components.

Each Component of the Whole Rice Plant	Length (mm)		Width (Outer Diameter) (mm)		Thickness (Inner Diameter) (mm)		Moisture Content (%)		Density (kg/m3)
	Average Value	Standard Deviation	Average Value	Standard Deviation	Average Value	Standard Deviation	Average Value	Standard Deviation	Minimum Value	Maximum Value	Average Value	Standard Deviation
Rice grains	5.86	0.16	2.10	0.08	1.86	0.07	22.67	3.22	1050	1330	1200	14
Main stem	837.62	10.41	5.18	0.46	4.02	0.41	66.54	7.67	182	224	216	12
Primary stem	76.73	4.22	/	/	0.98	0.02	52.65	5.33	179	184	181	9
Secondary stem	5.08	0.17	/	/	0.96	0.02
The whole plant	902.67	12.34	/	/	/	/	/	/	/	/	/	/

,

Table A3. Friction coefficient between various components of the whole rice plant and contact materials.

Parameter		Numerical Value
		Minimum Value	Maximum Value	Average Value	Standard Deviation
Static friction coefficient	Rice grain–rice grain	0.30	0.66	0.48	0.08
	Rice grain–main stem	0.62	0.71	0.65	0.07
	Rice grain–branch	0.55	0.62	0.60	0.08
	Rice grain–steel	0.26	0.58	0.42	0.03
	Main stem–main stem	0.36	0.76	0.55	0.03
	Main stem–branch	0.52	0.66	0.54	0.05
	Main stem–steel	0.32	0.68	0.46	0.02
	Branch–branch	0.68	0.85	0.70	0.04
	Branch–steel	0.42	0.65	0.55	0.07
Rolling friction coefficient	Rice grain–rice grain	0.010	0.020	0.015	0.002
	Rice grain–main stem	0.012	0.036	0.020	0.004
	Rice grain–branch	0.012	0.024	0.018	0.002
	Rice grain–steel	0.005	0.015	0.010	0.001
	Main stem–main stem	0.010	0.060	0.040	0.001
	Main stem–branch	0.028	0.040	0.032	0.003
	Main stem–steel	0.010	0.050	0.030	0.001
	Branch–branch	0.038	0.056	0.046	0.006
	Branch–steel	0.030	0.044	0.038	0.003

,

Table A4. Elastic modulus of the whole rice plant.

Each Component of the Whole Rice Plant	Elastic Modulus
	Minimum Value	Maximum Value	Average Value	Standard Deviation
Rice grain	3515	3875	3640	128
Main stem	314	335	320	7.6
Branch	120	250	150	8.4

and

Table A5. Collision restitution coefficient of the whole rice.

Parameters	Value
	Minimum Value	Maximum Value	Average Value	Standard Deviation
Rice grain–rice grain	0.15	0.50	0.42	0.07
Rice grain–main stem	0.45	0.52	0.44	0.05
Rice grain–branch	0.38	0.58	0.42	0.11
Rice grain–steel	0.43	0.57	0.56	0.06
Main stem–main stem	0.10	0.46	0.25	0.06
Main stem–branch	0.40	0.53	0.44	0.09
Main stem–steel	0.24	0.52	0.37	0.06
Branch–branch	0.25	0.38	0.32	0.02
Branch–steel	0.32	0.44	0.39	0.08

.

2.2. The Discrete Element Mechanical Model

Two kinds of models were applied in the discrete element modeling process of the whole rice plant. The surface of the rice grains during the harvest period was smooth and noncohesive. The Hertz–Mindlin (no slip) model was used to construct the contact between rice grains and rice grains and between rice grains and other parts. In other words, the changes in the rice grain parameters, such as force, displacement, and velocity, were all caused by the deformation due to collision [31]. In addition to other parts of rice grains, the whole rice plant showed flexibility. In the modeling process, the Hertz–Mindlin with bonding model was used to construct main stem–main stem, main stem–branch, branch–branch, and branch–grain contacts [32].

In addition, when the rice stands upright in the field during the harvest period, the particle contact model of the fixing device at the bottom of the rice is also set to the Hertz–Mindlin with bonding model so that strong contact adhesion can be formed between the particles and ensure the upright characteristics of rice during the harvest period. The relevant discrete element mechanical model is shown in Figure 4. The relevant theories are shown in Appendix A.

In the bonding bond setting process of Hertz–Mindlin with the bonding model, the required parameters mainly include normal stiffness coefficient, tangential stiffness coefficient, critical normal stress, critical tangential stress, and particle bonding radius.

The parameters of the bonding model were determined mainly by early compression and shear tests.

$$\{\begin{array}{l}{K}_n=\frac{4}{3}{\left(\frac{1-v_i^2}{E_i}+\frac{1-v_j^2}{E_j}\right)}^{-1}{\left(\frac{r_i+r_j}{r_i{r}_j}\right)}^{-\frac{1}{2}}\\ K_s=\left(\frac{1}{2}~\frac{2}{3}\right)K_n\end{array}$$

(13)

where K_n is the normal stiffness coefficient, N/m; K_s is the tangential stiffness coefficient, N/m; v_i is the Poisson’s ratios of particles i; v_j is the Poisson’s ratios of the particles j; E_i is the elastic modulus of the particle i, MPa; E_j is the elastic modulus of the particle j, MPa; r is the radius of the particle i, mm; r_j is the radius of the particle j, mm.

In addition, the particle bonding radius is generally 1.2 to 2 times the particle radius, where the particle radius is set according to the geometric characteristics of the main stem and branch of rice.

2.3. Establishment of Discrete Element Model

Based on a combined approach of solid modeling and API coordinate positioning, a discrete element model of the whole rice plant was established. Based on the previously determined geometric characteristics of the whole rice plant, the plant and its bottom fixing device were divided into main stem, primary stem, secondary stem, grain, and fixing device. The CATIA V5R20 software (Dassault Systèmes, Paris, France) was used to create a three-dimensional model and imported into ANSYS Workbench (version: 18.0, ANSYS Inc., Canonsburg, PA, USA) for meshing. The main stem was divided into cubic grids of 5 mm × 5 mm × 5 mm, the primary and secondary stems into cubic grids of 1 mm × 1 mm × 1 mm, the grains into rectangular grids of 2 mm × 2 mm × 5 mm, and the fixing device into cubic grids of 10 mm × 10 mm × 10 mm. Five mesh files were imported into Fluent, and the user-defined functions (UDF) file was used to solve the coordinates of each grid’s center point and generate coordinate files. The particles in the five coordinate files were named to ensure consistency with the particle template names in EDEM software (version: 2018, Altair Inc., Troy, MI, USA). Finally, the five coordinate files were merged into one and imported into EDEM API. Particle templates were created for each component in EDEM. Since the particle coordinates in the template are relative coordinates, the positions of each particle can be determined based on the coordinates in the API. The single-sphere filling was applied to the main stem, primary and secondary stems, and the fixing device, while multisphere filling was used for the rice grains, as shown in Figure 3 (the multisphere coordinates are shown in

Table A6. Multi-spherical coordinates (relative coordinates of rice grains) and radii.

Order	X Coordinate (mm)	Y Coordinate (mm)	Z Coordinate (mm)	Physical Radius (mm)
1	0	0	−2	0.5
2	0	0	−1	1
3	0	0	0	1.2
4	0	0	1	1
5	0	0	2	0.2

). The relative magnification factors of the main stem, primary stem, secondary stem, fixture, and rice grain are 2.63, 0.53, 0.53, 1.33, and 5.26, respectively, and a total of 43,214 particles were required to fill all five components.

The EDEM API was imported into the particle factory, ensuring that the coordinate file was in the same folder as the API. Simultaneously, the divided five-part mesh grid was imported into the EDEM entity, setting the type as Physical. This effectively restricts the relative coordinates of the particles generated in the particle factory. Automatically the computational domain was updated. The contact between the particles was set as Hertz–Mindlin with bonding in the Physics settings. The time step was set as 5%, with a particle generation time of 0.001 s and a bonding generation time of 0.0012 s among particles. The total simulation time was set as 0.002 s, and the grid size was set as 3R_min. The resulting bonding key model is shown in Figure 5c.

The total number of bonding keys formed collectively was 53,473, including 5004 bonding keys among main stems, 1368 bonding keys between main stems and primary stems, 72 bonding keys between main stems and fixed device, 22,752 bonding keys among primary stems, 3024 bonding keys between primary stems and secondary stems, 576 bonding keys between primary stems and rice grains, 8064 bonding keys among secondary stems, 2016 bonding keys between secondary stems and rice grains, and 10,597 bonding keys among fixed device.

2.4. Parameter Calibration Test

2.4.1. Parameter Calibration of Discrete Element Model for Rice Grain

The main contact materials include rice grains–rice grains, rice grains–steel, rice grains–main stems, and rice grains–branches. In addition, the density and elastic modulus of the rice grains needed to be set, and for each contact material, the static friction coefficient, coefficient of rolling friction, and restitution coefficient needed to be determined. The calibration tests required for these parameters are extensive. The proposed approach in this study was to conduct Plackett–Burman screening experiments [33], steepest ascent experiments, and Box–Behnken optimization experiments using rice grains–steel contact as an example. By doing so, the density and elastic modulus of rice grains, as well as the static friction coefficient, rolling friction coefficient, and restitution coefficient of rice grains–rice grains and rice grains–steel, can be determined. The known parameters are then used to conduct the steepest ascent experiments and Box–Behnken optimization experiments on various contact materials to determine the discrete element parameters of rice grains with different contact materials. DEM simulations were carried out working on an HPE ProLiant XL190r computer with Intel^® Xeon^® Gold 5120 CPU @ 2.20 GHz and 512 GB RAM.

The physical parameters of the tested rice grains and the contact parameters between rice grains and steel are encoded as −1 and +1 levels, representing the minimum and maximum values, respectively, as shown in

Table 1. Test parameter of Plackett–Burman (rice grains).

Simulation Parameters	Level
	−1	+1
Static friction coefficient (rice grains–rice grains) x1	0.30	0.66
Static friction coefficient (rice grains–steel) x2	0.26	0.58
Rolling friction coefficient (rice grains–rice grains) x3	0.01	0.02
Rolling friction coefficient (rice grains–steel) x4	0.005	0.015
Restitution coefficient (rice grains–rice grains) x5	0.15	0.50
Restitution coefficient (rice grains–steel) x6	0.43	0.57
Density of rice grains x7 (kg/m3)	1050	1330
Elastic modulus of rice grains x8 (MPa)	3515	3875

.

The natural repose angle of rice grains was calibrated through bench tests and virtual simulation experiments, as shown in Figure 6. A vertical cylindrical container made of Q235 steel with a diameter of 50 mm and a height of 150 mm was placed on a horizontal surface. The z-axis direction was chosen as the gravity direction, with a value set at −9.81 m/s². The “rainfall method” was used to allow the rice grains to naturally accumulate inside the cylinder, starting from an initial velocity of zero. The total number of rice grains was set to 4000. Once the rice grains settled completely inside the cylinder, the cylinder was lifted at a speed of 0.05 m/s until the rice grains naturally dispersed and formed a stable heap under the influence of gravity. The simulation was stopped when the rice grain heap reached a completely static and stable state. The Rayleigh time step was determined by the software based on parameters such as the grain radius and density, and it was ultimately set at 20% with a grid size of 3R_min. The total simulation duration for each group was set to 3 s. The duration of each realistic simulation is approximately 11 h.

After the test, in order to accurately measure the natural repose angle of rice grains, the three-dimensional population heap of rice grains was captured using cameras. MATLAB software (version: R2017b, MathWorks Company, Natick, MA, USA) was used for image noise reduction, grayscale processing, and binarization to extract the contour curve of the rice grain heap from one side. The size of the natural repose angle was determined by the ratio of vertical and horizontal pixel values on the curve. Each test was repeated three times, and the average result was taken.

2.4.2. Calibration of Discrete Element Model Parameters for Main Stems and Branches

The main stems and branches exhibited flexible characteristics. Therefore, in the calibration process of flexible main stems or branches, the calibration indicators included not only the natural rest angle but also the deflection.

To begin, a test was conducted to examine the contact between the rice main stem and steel. This involved Plackett–Burman screening tests, steepest ascent tests, and Box–Behnken optimization tests, each performed separately. During this phase, parameters such as density, the elastic modulus of the main stems, static friction coefficients between the main stem–main stem, main stem–steel, rolling friction coefficient, restitution coefficient, normal stiffness coefficient, tangential stiffness coefficient, critical normal stress, critical tangential stress, and particle bonding radius of the main stems were determined. These established parameters were then utilized as known quantities to conduct the steepest ascent tests and Box–Behnken optimization tests on various contact materials, thereby determining the discrete element parameters between the main stems and different contact materials.

The physical parameters of the tested main stems, the contact parameters between the main stems and steel, and the bonding parameters were encoded as −1 and +1 levels, representing the minimum and maximum values, respectively, as shown in

Table 2. Test parameter of Plackett–Burman (main stems).

Simulation Parameters	Level
	−1	+1
Static friction coefficient (main stems—main stems) X1	0.36	0.76
Static friction coefficient (main stems—steel) X2	0.32	0.68
Coefficient of rolling friction (main stems—main stems) X3	0.01	0.06
Coefficient of rolling friction (main stems—steel) X4	0.01	0.05
Restitution coefficient (main stems—main stems) X5	0.10	0.46
Restitution coefficient (main stems—steel) X6	0.24	0.52
Density of main stems X7 (kg/m3)	182	224
Elastic modulus of main stems X8 (MPa)	314	335
Normal stiffness coefficient X9 (N/m)	1.34 × 1012	1.60 × 1012
Tangential stiffness coefficient X10 (N/m)	1.22 × 1012	1.50 × 1012
Critical normal stress X11 (Pa)	1.15 × 1011	1.30 × 1011
Critical tangential stress X12 (Pa)	1.20 × 1011	1.32 × 1011
Bonding radius X13 (mm)	3	5

.

The natural repose angle and deflection of the main stem were calibrated through bench tests and virtual simulation tests, as shown in Figure 7. Each test was repeated three times, and the average value was taken as the result.

• Natural rest angle calibration: A main stem model with a pitch of 20 mm was established. A cylindrical container made of Q235 steel material with the following physical properties: density of 7850 kg/m³, elastic modulus of 212 MPa, and Poisson’s ratio of 0.31 was placed vertically on a horizontal plane. The diameter of the cylindrical container was 50 mm, and its height was 150 mm. The z-axis direction was selected as the gravity direction, with a value set to −9.81 m/s². The “rainfall method” was employed to allow the main stems to naturally pile up inside the cylindrical container, starting from zero initial velocity. A total of 300 main stems were used in the simulation. Once the main stems had completely settled inside the cylindrical container, the container was lifted at a velocity of 0.05 m/s until the main stems naturally scattered and formed a stable pile under the influence of gravity. The simulation was stopped when the pile of main stems reached a completely static and stable state. The Rayleigh time step was set to 20%, and the grid size was 3R_min. Each simulation group had a total duration of 3 s. The duration of each realistic simulation is approximately 5 h.
• Deflection calibration: The deflection of rice main stems was calibrated using a three-point bending test. In bench tests and simulation tests, it is necessary to ensure the uniformity of the internodes of the stem and the mass and material of the calibrated object. A main stem model with a pitch of 130 mm was established and placed in the middle of a support model. A wedge-shaped steel block weighing 50 g was placed in the middle of the main stem. The Rayleigh time step was set to 20%, and the grid size was 3R_min. Each simulation group had a total duration of 3 s. During the compression process, the deflection of the main stem was monitored. The duration of each realistic simulation is approximately 20 min.

3. Results and Discussion

3.1. The Calibration of Discrete Element Parameters for Rice Grains

The Plackett–Burman test plan and results for the calibration of discrete element parameters of rice grains are shown in

Table 3. The Plackett–Burman test plan and results.

Order	Test Factors								Index
	x1	x2	x3	x4	x5	x6	x7	x8	Natural Repose Angle y1/(°)
1	1	1	−1	1	1	1	−1	−1	20.65
2	−1	1	1	−1	1	1	1	−1	21.32
3	1	−1	1	1	−1	1	1	1	23.90
4	−1	1	−1	1	1	−1	1	1	23.54
5	−1	−1	1	−1	1	1	−1	1	23.77
6	−1	−1	−1	1	−1	1	1	−1	25.26
7	1	−1	−1	−1	1	−1	1	1	23.65
8	1	1	−1	−1	−1	1	−1	1	22.70
9	1	1	1	−1	−1	−1	1	−1	20.85
10	−1	1	1	1	−1	−1	−1	1	23.54
11	1	−1	1	1	1	−1	−1	−1	22.93
12	−1	−1	−1	−1	−1	−1	−1	−1	24.47

.

The significance analysis of the regression model for the natural rest angle y₁ is conducted, and the results are shown in

Table 4. Significance analysis of Plackett–Burman test results.

Source	Sum of Squares	Degree of Freedom	Mean Square	F-Value	p-Value	Significance
Model	22.00	8	2.75	11.59	0.0344	*
x1	4.34	1	4.34	18.32	0.0234	*
x2	10.79	1	10.79	45.50	0.0067	**
x3	1.31	1	1.31	5.51	0.1006
x4	0.78	1	0.78	3.29	0.1674
x5	1.97	1	1.97	8.30	0.0635
x6	0.16	1	0.16	0.67	0.4733
x7	0.018	1	0.018	0.074	0.8028
x8	2.63	1	2.63	11.10	0.0447	*
Residual	0.71	3	0.24
Cor Total	22.71	11

** Represents highly significant (p < 0.01); * represents significant (p < 0.05).

. The regression model for y₁ is as follows:

$$y_1=23.05-0.60x_1-0.95x_2-0.33x_3+0.26x_4-0.40x_5-0.11x_6+0.038x_7+0.47x_8$$

(14)

The model had a p-value of less than 0.05, indicating the significant influence of the main effects model. The fitted regression equation aligned well with the actual situation. Among them, x₂ had an highly significant impact on the natural repose angle y₁, while x₁ and x₈ had a significant impact. The other factors were not significant. Based on the regression equation and significance analysis, the primary and secondary factors influencing the natural repose angle y₁ were x₂, x₁, x₈, x₅, x₃, x₄, x₆, and x₇. To determine their effects on the index, the steepest ascent test was conducted on x₁, x₂, and x₈, which significantly affected the natural repose angle.

The natural rest angle y₁ and the relative error of the angle of repose $\eta$ were chosen as the test and evaluation indicators. The steepest ascent test plan and results are shown in

Table 5. The steepest ascent test plan and results.

Order	Factors			Index
	Static Friction Coefficient (Rice Grains–Rice Grains) x1	Static Friction Coefficient (Rice Grains–Steel) x2	Elastic Modulus of Rice Grains x8 (MPa)	Natural Repose Angle y1 (°)	Relative Error of the Natural Repose Angle $\mathit{\eta }$ (%)
1	0.30	0.26	3515	20.34	13.37
2	0.39	0.34	3605	22.40	4.60
3	0.48	0.42	3695	23.17	1.32
4	0.57	0.50	3785	25.46	8.43
5	0.66	0.58	3875	27.23	15.97

. The formula for calculating the relative error of the repose angle $\eta$ is as follows:

$$\eta =\frac{\left|y-y_1\right|}{y}\times 100\%$$

(15)

where $\eta$ is the relative error of the angle of repose, %; y₁ is the natural repose angle during simulation, (°); and y is the natural repose angle during test, (°).

As the static friction coefficient (rice grains–rice grains) x₁, the static friction coefficient (rice grains–steel) x₂, and the elastic modulus of rice grains x₈ gradually increased, the natural repose angle y₁ also increased gradually, while the relative error of the angle of repose initially decreased and then increased. When the values of these three factors were set to trial order 3 in the steepest ascent test, the relative error of the angle of repose was minimized. Therefore, for the Box–Behnken test, the parameters of trial order 3 from the steepest ascent test were chosen as the middle level (0), while trial orders 2 and 4 were set as the low (−1) and high (+1) levels, respectively, to investigate the effects of these three test factors on the natural repose angle. The levels encoding table for the test factors are shown in

Table 6. Levels encoding table for the Box–Behnken test factors.

Codes	Test Factors
	Static Friction Coefficient (Rice Grains–Rice Grains) x1	Static Friction Coefficient (Rice Grains–Steel) x2	Elastic Modulus of Rice Grains x8 (MPa)
−1	0.39	0.34	3605
0	0.48	0.42	3695
1	0.57	0.50	3785

, and the plan and results of the Box–Behnken test are presented in

Table 7. Box–Behnken test plan and results.

Order	Test Factors			Index
	Static Friction Coefficient (Rice Grains–Rice Grains) x1	Static Friction Coefficient (Rice Grains–Steel) x2	Elastic Modulus of Rice Grains x8 (MPa)	Natural Repose Angle y1 (°)
1	−1	−1	0	23.65
2	1	−1	0	24.17
3	−1	1	0	22.89
4	1	1	0	26.54
5	−1	0	−1	21.03
6	1	0	−1	22.78
7	−1	0	1	23.97
8	1	0	1	24.12
9	0	−1	−1	23.05
10	0	1	−1	24.66
11	0	−1	1	25.7
12	0	1	1	26.03
13	0	0	0	24.08
14	0	0	0	24.36
15	0	0	0	24.11
16	0	0	0	24.07
17	0	0	0	23.95

.

The test results were analyzed using Design-Expert 8.0.6 software (Stat-Ease, Minneapolis, MN, USA) to investigate the effects of the three test factors on the natural rest angle. The analysis of the variance table for the test factors is shown in

Table 8. Analysis of variance for the factors on the test indicator.

Source	Sum of Squares	Degree of Freedom	Mean Square	F-Value	p-Value	Significance
Model	25.83	9	2.87	26.55	0.0001	**
x1	4.61	1	4.61	42.61	0.0003	**
x2	1.58	1	1.58	14.57	0.0066	**
x8	8.61	1	8.61	79.66	<0.0001	**
x1x2	2.45	1	2.45	22.66	0.0021	**
x1x8	0.64	1	0.64	5.92	0.0452	*
x2x8	0.41	1	0.41	3.79	0.0926
$x_1^2$	2.99	1	2.99	27.70	0.0012	**
$x_2^2$	4.57	1	4.57	42.27	0.0003	**
$x_8^2$	0.37	1	0.37	3.41	0.1074
Residual	0.76	7	0.11
Cor Total	26.59	16

** Represents highly significant (p < 0.01); * represents significant (p < 0.05).

.

The analysis of variance for the factors on the natural repose angle reveals that factors x₁, x₂, x₈, x₁x₂, $x_1^2$, and $x_2^2$ exhibit a significance level of less than 0.01, indicating their highly significant impact. Factor x₁x₈ demonstrates a significance level of less than 0.05, signifying its significant impact. However, factors x₂x₈ and $x_8^2$ exhibited a significance level greater than 0.05, suggesting they were insignificance.

The regression model for the natural repose angle is highly significant, indicating a good fit. Using the natural repose angle y₁ as the response function and the factors as independent variables, the fitted regression equation is as follows:

$$y_1=24.11+0.76x_1+0.44x_2+1.04x_8+0.78x_1{x}_2-0.4x_1{x}_8-0.32x_2{x}_8-0.84x_1^2+1.04x_2^2-0.3x_8^2$$

(16)

On this basis, the Design-Expert 8.0.6 software was utilized to establish response surface plots depicting the effects of the static friction coefficient (rice grains–rice grains), static friction coefficient (rice grains–steel), and elastic modulus of rice grains on the natural repose angle. The response surface plot is shown in Figure 8.

When the static friction coefficient (rice grains–steel) was held constant, the natural repose angle gradually increased with an increasing static friction coefficient (rice grains–rice grains), and also increased gradually with an increasing elastic modulus of rice grains. When the static friction coefficient (rice grains–rice grains) remained constant, the natural repose angle initially decreased and then increased with an increasing static friction coefficient (rice grains–steel), while it gradually increased with an increasing elastic modulus of rice grains. When the elastic modulus of rice grains was kept constant, the natural repose angle gradually increased with an increasing static friction coefficient (rice grains–rice grains) and initially decreased and then increased with an increasing static friction coefficient (rice grains–steel).

Based on the above analysis, the significance of the factors affecting the natural repose angle can be ranked as follows: elastic modulus > static friction coefficient (rice grains–rice grains) > static friction coefficient (rice grains–steel).

Using the optimization module in Design Expert 8.0.6 software, the parameters of the three factors were optimized based on the target average value of the natural repose angle, which was 23.48°. The optimized values obtained were a static friction coefficient (rice grains–rice grains) of 0.39, a static friction coefficient (rice grains–steel) of 0.50, and an elastic modulus of rice grains of 3707 MPa.

With the optimized values of the significant parameters and the average values of the nonsignificant parameters, three validation experiments were conducted, resulting in natural angles of repose of 22.96°, 23.23°, and 23.78°, with an average of 23.32°. The relative error compared to the actual measured values was 0.68%, validating the accuracy of the related discrete element simulation parameter settings.

The determined parameters for the contact between rice grains and different materials were used as known quantities to calibrate the static friction coefficient, rolling friction coefficient, and restitution coefficient between rice grains–main stems and rice grains–branches. Since there are only three unknown parameters for each contact material, the Plackett–Burman screening test method was not used to identify significant influencing factors. Instead, the steepest ascent test and Box–Behnken optimization test were directly conducted to narrow down the range of values for the three unknown factors and find the optimal combination of parameters.

The test methods for determining the contact parameters between rice grains and different materials were the same, and the detailed analysis process is not repeated here. The final obtained discrete element parameters for the contact between rice grains and different materials are presented in

Table 9. Numerical values of discrete element parameters for the contact between rice grains and different materials.

Parameters	Numerical Value
Density of rice grains (kg/m3)	1200
Elastic modulus of rice grains (MPa)	3707
Poisson’s ratio of rice grains *	0.25
Static friction coefficient (rice grains–rice grains)	0.39
Static friction coefficient (rice grains–steel)	0.50
Static friction coefficient (rice grains–main stems)	0.65
Static friction coefficient (rice grains–branches)	0.60
Rolling friction coefficient (rice grains–rice grains)	0.015
Rolling friction coefficient (rice grains–steel)	0.010
Rolling friction coefficient (rice grains–main stems)	0.020
Rolling friction coefficient (rice grains–branches)	0.020
Restitution coefficient (rice grains–rice grains)	0.42
Restitution coefficient (rice grains–steel)	0.56
Restitution coefficient (rice grains–main stems)	0.395
Restitution coefficient (rice grains–branches)	0.420

* Denotes values referenced from the relevant literature [30].

.

3.2. Calibration of Discrete Element Parameters for Stems

The Plackett–Burman test plan and results for the calibration of discrete element parameters for stems are presented in

Table 10. Plackett–Burman test plan and results.

Order	Test Factors													Indexes
	X1	X2	X3	X4	X5	X6	X7	X8	X9	X10	X11	X12	X13	Natural Repose Angle Y1 (°)	Deflection Y2 (mm)
1	1	1	−1	−1	1	1	1	1	−1	1	−1	1	−1	38.98	40.93
2	−1	1	1	−1	−1	1	1	1	1	−1	1	−1	1	35.34	43.40
3	1	−1	1	1	−1	−1	1	1	1	1	−1	1	−1	35.23	42.00
4	1	1	−1	1	1	−1	−1	1	1	1	1	−1	1	31.53	43.14
5	−1	1	1	−1	1	1	−1	−1	1	1	1	1	−1	35.92	41.96
6	−1	−1	1	1	−1	1	1	−1	−1	1	1	1	1	33.95	42.80
7	−1	−1	−1	1	1	−1	1	1	−1	−1	1	1	1	37.3	44.00
8	−1	−1	−1	−1	1	1	−1	1	1	−1	−1	1	1	32.1	44.62
9	1	−1	−1	−1	−1	1	1	−1	1	1	−1	−1	1	34.55	42.97
10	−1	1	−1	−1	−1	−1	1	1	−1	1	1	−1	−1	32.07	40.67
11	1	−1	1	−1	−1	−1	−1	1	1	−1	1	1	−1	31.85	43.44
12	−1	1	−1	1	−1	−1	−1	−1	1	1	−1	1	1	31.12	45.13
13	1	−1	1	−1	1	−1	−1	−1	−1	1	1	−1	1	32.64	42.90
14	1	1	−1	1	−1	1	−1	−1	−1	−1	1	1	−1	36.56	41.78
15	1	1	1	−1	1	−1	1	−1	−1	−1	−1	1	1	39.5	42.35
16	1	1	1	1	−1	1	−1	1	−1	−1	−1	−1	1	33.8	43.78
17	−1	1	1	1	1	−1	1	−1	1	−1	−1	−1	−1	34.76	40.25
18	−1	−1	1	1	1	1	−1	1	−1	1	−1	−1	−1	31.67	39.70
19	1	−1	−1	1	1	1	1	−1	1	−1	1	−1	−1	38.52	41.94
20	−1	−1	−1	−1	−1	−1	−1	−1	−1	−1	−1	−1	−1	32.12	40.66

.

The significance analysis of the regression models for natural repose angle Y₁ and deflection Y₂ is presented in

Table 11. Analysis of variance for the Plackett–Burman test result.

Index	Source	Sum of Squares	Degree of Freedom	Mean Square	F-Value	p-Value	Significance
Y1	Model	120.35	13	9.26	4.23	0.0434	*
	X1	14.13	1	14.13	6.45	0.0441	*
	X2	4.66	1	4.66	2.13	0.1951
	X3	0	1	0	0	0.9780
	X4	0.02	1	0.02	0	0.9273
	X5	13.33	1	13.33	6.09	0.0486	*
	X6	8.80	1	8.80	4.02	0.0918
	X7	47.71	1	47.71	21.79	0.0034	**
	X8	4.77	1	4.77	2.18	0.1903
	X9	2.94	1	2.94	1.34	0.2905
	X10	10.07	1	10.07	4.60	0.0757
	X11	0.17	1	0.17	0.078	0.7892
	X12	12.03	1	12.03	5.49	0.0576
	X13	1.71	1	1.71	0.78	0.4107
	Residual	13.14	6	2.19
	Cor Total	133.49	19
Y2	Model	38.25	13	2.94	6.03	0.0184	*
	X1	0.22	1	0.22	0.45	0.5264
	X2	0.12	1	0.12	0.26	0.6310
	X3	0.55	1	0.55	1.13	0.3287
	X4	0.023	1	0.023	0.047	0.8349
	X5	1.20	1	1.20	2.46	0.1677
	X6	0.018	1	0.018	0.037	0.8540
	X7	1.72	1	1.72	3.52	0.1097
	X8	0.41	1	0.41	0.85	0.3920
	X9	4.25	1	4.25	8.72	0.0255	*
	X10	0.83	1	0.83	1.71	0.2393
	X11	0.68	1	0.68	1.40	0.2809
	X12	4.67	1	4.67	9.57	0.0213	*
	X13	23.54	1	23.54	48.28	0.0004	**
	Residual	2.93	6	0.49
	Cor Total	41.17	19

** Represents highly significant (p < 0.01); * represents significant (p < 0.05).

. The regression model for Y₁ is as follows:

$$\begin{array}{l}{Y}_1=34.48+0.84X_1+0.48X_2-0.032X_4+0.82X_5+0.66X_6+1.54X_7-0.49X_8-0.38X_9\\ \begin{array}{cc}& \end{array}-0.71X_{10}+0.092X_{11}+0.78X_{12}-0.29X_{13}\end{array}$$

(17)

The regression model for Y₁ is as follows:

$$\begin{array}{l}{Y}_2=42.42+0.11X_1-0.079X_2-0.17X_3+0.034X_4-0.24X_5-0.03X_6-0.29X_7+0.14X_8\\ \begin{array}{cc}& \end{array}+0.46X_9-0.20X_{10}+0.18X_{11}+0.48X_{12}+1.08X_{13}\end{array}$$

(18)

According to the significance analysis, the regression model for the indicator of natural repose angle Y₁ had a p-value less than 0.05, indicating that the main effect model was significant and the fitted regression equation aligns well with the actual situation. Among them, X₇ had an highly significant impact on the natural repose angle Y₁, while X₁ and X₅ had a significant impact. The other factors were not significant. Based on the regression equation and significance analysis, the primary and secondary factors influencing the natural repose angle Y₁ were X₁, X₅, X₁₂, X₁₀, X₆, X₈, X₂, X₉, X₁₃, X₁₁, X₄, and X₃. To determine their effects on the index, the steepest ascent test was conducted on X₁, X₅, and X₇, which significantly affected the natural repose angle.

According to the significance analysis, the regression model for the indicator of deflection Y₂ has a p-value less than 0.05, indicating that the main effect model is significant and the fitted regression equation aligns well with the actual situation. Among them, X₁₃ had an highly significant impact on the deflection Y₂, while X₉ and X₁₂ had a significant impact. The other factors were not significant. Based on the regression equation and significance analysis, the primary and secondary factors influencing the deflection Y₂ were X₁₃, X₁₂, X₉, X₇, X₅, X₁₀, X₁₁, X₃, X₈, X₁, X₂, X₄, and X₆. To determine their effects on the index, the steepest ascent test was conducted on X₉, X₁₂, and X₁₃, which significantly affected the deflection.

The static friction coefficient (main stems–main stems) X₁, restitution coefficient (main stems–main stems) X₅, and stem density X₇ have a significant impact on the natural repose angle Y₁. The normal stiffness coefficient X₉, critical shear stress X₁₂, and bonding radius X₁₃ have a significant influence on the deflection Y₂. Therefore, the steepest ascent test was conducted for each of these six factors. The natural repose angle Y₁ and the relative error of the angle of repose $\eta$ were used as the test and evaluation index for the calibration of the natural repose angle. The deflection Y₂ and the relative error of the deflection $\delta$ were used as the test and evaluation index for the calibration of the deflection. The formula for calculating the relative error of the deflection $\delta$ is as follows:

$$\delta =\frac{\left|Y_2-Y\right|}{Y_2}\times 100\%$$

(19)

where $\delta$ is the relative error of the deflection, %; Y₂ is simulated deflection, mm; and Y is tested deflection, mm.

The test plan and results are shown in

Table 12. Steepest ascent test plan and results of natural repose angle.

Order	Test Factors			Indexes
	Static Friction Coefficient (Main Stems–Main Stems) X1	Restitution Coefficient (Main Stems–Main Stems) X5	Density of Main Stems X7 (kg/m3)	Natural Repose Angle Y1 (°)	The Relative Error of the Natural Repose Angle $\mathit{\eta }$ (%)
1	0.36	0.10	182	28.15	18.21
2	0.46	0.19	192.5	29.36	14.70
3	0.56	0.28	203	32.17	6.54
4	0.66	0.37	213.5	35.16	2.15
5	0.76	0.46	224	38.40	11.56

and

Table 13. The steepest ascent test plan and results of deflection.

Order	Factors			Indexes
	Normal Stiffness Coefficient X9 (N/m)	Critical Tangential Stress X12 (Pa)	Bonding Radius X13 (mm)	Deflection Y2 (mm)	The Relative Error of the Deflection $\mathit{\delta }$ (%)
1	1.34 × 1012	1.20 × 1011	3	45.62	7.58
2	1.405 × 1012	1.23 × 1011	3.5	43.27	2.01
3	1.47 × 1012	1.26 × 1011	4	39.62	6.56
4	1.535 × 1012	1.29 × 1011	4.5	37.24	12.17
5	1.60 × 1012	1.32 × 1011	5	35.45	16.39

.

As the static coefficient of friction (main stems–main stems) X₁, restitution coefficient (main stems–main stems) X₅, and main stem density X₇ gradually increase, the natural repose angle Y₁ also increases gradually, while the relative error of the angle of repose $\eta$ initially decreases and then increases. When the values of these three factors were set to order 4 in the test plan, the relative error of the angle of repose was minimized. Therefore, the parameters from the steepest ascent test corresponding to order 4 were taken as the middle level 0, and the factors for orders 3 and 5 were set as low −1 and high +1, respectively, for the Box–Behnken test. The level encoding table for the Box–Behnken test is shown in

Table 14. Level encoding table for the Box–Behnken test of the natural repose angle.

Codes	Test Factors
	Static Friction Coefficient (Main Stems–Main Stems) X1	Restitution Coefficient (Main Stems–Main Stems) X5	Density of Main Stem X7 (kg/m3)
−1	0.56	0.28	203
0	0.66	0.37	213.5
1	0.76	0.46	224

, and the Box–Behnken test plan and results are presented in

Table 15. Box–Behnken test plan and results of natural repose angle.

Order	Test Factors			Index
	Static Friction Coefficient (Main Stems–Main Stems) X1	Restitution Coefficient (Main Stems–Main Stems) X5	Density of Main Stems X7 (kg/m3)	Natural Repose Angle Y1 (°)
1	−1	−1	0	36.42
2	1	−1	0	35.60
3	−1	1	0	32.70
4	1	1	0	35.42
5	−1	0	−1	34.45
6	1	0	−1	37.58
7	−1	0	1	35.61
8	1	0	1	33.76
9	0	−1	−1	39.52
10	0	1	−1	38.47
11	0	−1	1	38.62
12	0	1	1	35.03
13	0	0	0	36.56
14	0	0	0	35.98
15	0	0	0	36.12
16	0	0	0	35.26
17	0	0	0	35.44

.

As the normal stiffness coefficient X₉, critical tangential stress X₁₂, and bonding radius X₁₃ gradually increase, the deflection Y₂ decreases gradually, while the relative error of the deflection $\delta$ initially decreases and then increases. When the values of these three factors were set to order 3 in the test plan, the relative error of the deflection was minimized. Therefore, the factors from the steepest ascent test corresponding to order 3 were taken as the middle level 0, and the factors for orders 2 and 4 were set as low −1 and high +1, respectively, for the Box–Behnken test. The level encoding table for the Box–Behnken test is shown in

Table 16. Level encoding table for the Box–Behnken test of the deflection.

Codes	Test Factor
	Normal Stiffness Coefficient X9 (N/m)	Critical Tangential Stress X12 (Pa)	Bonding Radius X13 (mm)
−1	1.405 × 1012	1.23 × 1011	3.5
0	1.47 × 1012	1.26 × 1011	4
1	1.535 × 1012	1.29 × 1011	4.5

, and the Box–Behnken test plan and results are shown in

Table 17. Box–Behnken test plan and results of the deflection.

Order	Test Factors			Index
	Normal Stiffness Coefficient X9 (N/m)	Critical Tangential Stress X12 (Pa)	Bonding Radius X13 (mm)	Deflection Y2 (mm)
1	−1	−1	0	48.45
2	1	−1	0	43.25
3	−1	1	0	44
4	1	1	0	43.25
5	−1	0	−1	47.57
6	1	0	−1	41.02
7	−1	0	1	37.93
8	1	0	1	42.95
9	0	−1	−1	45.17
10	0	1	−1	44.2
11	0	−1	1	46.1
12	0	1	1	38.48
13	0	0	0	46.84
14	0	0	0	44.82
15	0	0	0	45.5
16	0	0	0	45.61
17	0	0	0	46.45

.

The test results were analyzed using Design-Expert 8.0.6 software to investigate the influence of test factors on the natural repose angle and deflection. The analysis of variance for the test factors is shown in

Table 18. Analysis of variance for the test factors on the test indexes.

Index	Source	Sum of Squares	Degree of Freedom	Mean Square	F-Value	p-Value	Significance
$Y_1$	Model	47.52	9	5.28	23.44	0.0002	**
	X1	1.26	1	1.26	5.61	0.0497	*
	X5	9.12	1	9.12	40.47	0.0004	**
	X7	6.13	1	6.13	27.19	0.0012	**
	X1X5	3.13	1	3.13	13.91	0.0074	**
	X5X7	6.20	1	6.20	27.52	0.0012	**
	X1X7	1.61	1	1.61	7.16	0.0317	*
	X12	12.15	1	12.15	53.92	0.0002	**
	X52	3.12	1	3.12	13.87	0.0074	**
	X72	5.83	1	5.83	25.87	0.0014	**
	Residual	1.58	7	0.23
	Cor Total	49.10	16
$Y_2$	Model	127.38	9	14.15	12.04	0.0017	**
	X9	6.99	1	6.99	5.95	0.0448	*
	X10	21.26	1	21.26	18.08	0.0038	**
	X13	19.53	1	19.53	16.61	0.0047	**
	X9X10	4.95	1	4.95	4.21	0.0793
	X9X13	33.47	1	33.47	28.46	0.0011	**
	X10X13	11.06	1	11.06	9.40	0.0182	*
	X92	5.22	1	5.22	4.44	0.0732
	X102	0	1	0	0	0.9902
	X132	23.52	1	23.52	20.00	0.0029	**
	Residual	8.23	7	1.18
	Cor Total	135.61	16

** Represents highly significant (p < 0.01); * represents significant (p < 0.05).

.

The analysis of variance for the factors on the natural repose angle indicated that factors X₅, X₇, X₁X₅, X₅X₇, X₁², X₅², and X₇² had a significance level of less than 0.01, indicating their highly significant impact. Factor X₁ and X₁X₇ had a significance level of less than 0.05, signifying their significant influence.

The regression model for the natural repose angle is highly significant, indicating a good fit. Using the natural repose angle Y₁ as the response function and the factors as independent variables, the fitted regression equation is as follows:

$$\begin{array}{l}{Y}_1=35.87+0.4X_1-1.07X_5-0.88X_7+0.88X_1{X}_5-1.24X_1{X}_7-0.63X_5{X}_7\\ \begin{array}{cc}& \end{array}-1.7X_1^2+0.86X_2^2+1.18X_3^2\end{array}$$

(20)

The analysis of variance for the factors on the deflection indicated that factors X₁₀, X₁₃, X₉X₁₃, X₁₃², and X₇² had a significance level of less than 0.01, indicating their highly significant impact. Factor X₉ and X₁₀X₁₃ had a significance level of less than 0.05, signifying their significant influence.

The regression model for the deflection is highly significant, indicating a good fit. Using the natural repose angle Y₂ as the response function and the factors as independent variables, the fitted regression equation is as follows:

$$\begin{array}{l}{Y}_2=45.84-0.94X_9-1.63X_{10}-1.56X_{13}+1.11X_9{X}_{10}+2.89X_9{X}_{13}\\ \begin{array}{cc}& \end{array}-1.66X_{10}{X}_{13}-1.11X_9^2-2.36X_{13}^2\end{array}$$

(21)

Based on this, Design-Expert 8.0.6 software was used to create response surface plots for the impact of the static friction coefficient (main stems–main stems), restitution coefficient (main stems–main stems), and density of main stems on the natural repose angle. Additionally, response surface plots were generated for the influence of the normal stiffness coefficient, critical tangential stress, and bonding radius on the deflection, as shown in Figure 9.

When the restitution coefficient (main stems–main stems) was held constant, the natural repose angle initially increased and then decreased with an increasing static friction coefficient (main stems–main stems), and it initially decreased and then increased with an increasing density of main stems. When the static friction coefficient (main stems–main stems) remained constant, the natural repose angle gradually decreased with an increasing restitution coefficient (main stems–main stems), and it initially decreased and then increased with an increasing main stem density. When the density of the main stems was kept constant, the natural repose angle gradually increased with an increasing static friction coefficient (main stems–main stems), and it initially decreased and then increased with an increasing restitution coefficient (main stems–main stems).

When the critical tangential stress was held constant, the deflection gradually decreased with an increasing normal stiffness coefficient and an increasing bonding radius. When the normal stiffness coefficient remained constant, the deflection gradually decreased with an increasing critical tangential stress, while it initially increased and then decreased with an increasing bonding radius. When the bonding radius was kept constant, the deflection gradually decreased with an increasing normal stiffness coefficient and a decreasing critical tangential stress.

Using the optimization module in Design Expert 8.0.6 software, the relevant parameter values for the three factors were optimized with the target of achieving an average natural repose angle of 34.42°. The optimized values obtained were a static friction coefficient (main stems–main stems) of 0.73, a restitution coefficient (main stems–main stems) of 0.415 (rounded to 0.42), and a main stem density of 222 kg/m³. Taking these optimized values for the significant parameters, the other insignificant parameters were set to their average values for three validation tests. The resulting natural angles of repose were 33.68°, 32.71°, and 34.02°, with an average of 33.47°. The relative error compared to the actual measured values was 2.76%.

Using an average deflection value of 42.40 mm as the target, the relevant parameter values for the three factors were optimized. The optimized values obtained were a normal stiffness coefficient of 1.516 × 10¹² N/m, a critical tangential stress of 1.275 × 10¹¹ Pa, and a bonding radius of 4.4 mm. Taking these optimized values for the significant parameters, the other insignificant parameters were set to their average values for three validation tests. The resulting deflections were 43.18 mm, 42.76 mm, and 44.23 mm, with an average of 43.39 mm. The relative error compared to the actual measured values was 2.33%. This study validates the accuracy of the related discrete element simulation parameter settings.

The discrete element parameters for the interaction between the main stems and branches were calibrated using the determined parameters for the main stems–branches contact. Since there are only three unknown parameters for a specific contact material, the Plackett–Burman experimental method was not used for screening significant influencing factors. Instead, the steepest ascent test and Box–Behnken optimization test were directly conducted to narrow down the range of values for the three unknown factors and determine the optimal combination of parameters. The experimental methods for determining the contact parameters between the main stem and branches for different materials were the same, and the detailed analysis process is not reiterated here. The resulting discrete element parameters for the main stems, branches, and different materials are presented in

Table 19. Numerical values of discrete element parameters for the interaction between the main stems, branches, and different materials.

Parameter	Numerical Value
Density of main stems (kg/m3)	222
Elastic modulus of main stems (MPa)	320
Poisson’s ratio of main stems *	0.40
Static friction coefficient (main stems–main stems)	0.73
Static friction coefficient (main stems–steel)	0.46
Static friction coefficient (main stems–branches)	0.56
Rolling friction coefficient (main stems–main stems)	0.040
Rolling friction coefficient (main stems–steel)	0.030
Rolling friction coefficient (main stems–branches)	0.032
Restitution coefficient (main stems–main stems)	0.42
Restitution coefficient (main stems–steel)	0.47
Restitution coefficient (main stems–branches)	0.44
Normal stiffness coefficient of the main stems (N/m)	1.516 × 1012
Tangential stiffness coefficient of the main stems (N/m)	1.516 × 1012
Critical normal stress of the main stems (Pa)	1.275 × 1011
Critical tangential stress of the main stems (Pa)	1.275 × 1011
Bonding radius of the main stems (mm)	4.4
Density of branches (kg/m3)	181
Elastic modulus of branches (MPa)	150
Poisson’s ratio of branches *	0.40
Static friction coefficient (branches–branches)	0.75
Static friction coefficient (branches–steel)	0.55
Rolling friction coefficient (branches–branches)	0.046
Rolling friction coefficient (branches–steel)	0.035
Restitution coefficient (branches–branches)	0.32
Restitution coefficient (branches–steel)	0.36
Normal stiffness coefficient of the branches (N/m)	1.5 × 1010
Tangential stiffness coefficient of the branches (N/m)	1.0 × 1010
Critical normal stress of the branches (Pa)	5.0 × 108
Critical tangential stress of the branches (Pa)	5.0 × 108
Bonding radius of the branches (mm)	0.7

* Denotes numerical values referenced from the relevant literature [34].

.

3.3. Discussion

This study selected the widely cultivated “Longjing 29” rice variety from the Northeast region of China as the sample. It was divided into main stems, primary stems, secondary stems, and grains to investigate their respective physical parameters. The typical geometric characteristics, moisture content, and density were determined. The static friction coefficient and rolling friction coefficient were measured using a self-built friction coefficient testing platform and high-speed camera technology. The elastic modulus of the whole rice plant was determined through compression and tangential tests. The collision coefficient of multiple rice grains was accurately measured using the impact collision method. Plackett–Burman screening tests, steepest ascent tests, and Box–Behnken optimization tests were conducted to optimize the parameters, establishing accurate discrete element model parameters for each component of the whole rice plant.

In the process of measuring the rolling friction coefficient, certain errors may occur because of the following reasons:

• Rice grains, main stems, and branches have irregular shapes, leading to a combination of sliding and rolling on the material surface. It is difficult to distinguish the specific motion state during high-speed imaging observation.
• Small bouncing of rice components may occur on the grain board or branch board.
• Partial detachment of the material from the material surface can occur during the rolling process. These conditions can result in inaccuracies in the measurement of the rolling friction coefficient, necessitating multiple repetitions of the experiment.

At present, there is no established standard method or test device for measuring the rolling friction coefficient in agricultural materials, which leads to differences when applying different methods to measure the rolling friction coefficient for the same material. In future research, it is important to focus on classifying different agricultural materials, such as small-sized regular materials, large-sized irregular materials, and low-density hollow materials, in order to explore the optimal method for determining the dynamic coefficient of friction specific to each material. Additionally, auxiliary devices can be utilized to load the material properties without altering the contact characteristics, enabling effective measurement of the rolling friction coefficient in agricultural materials.

Because of the small size and light weight of rice grains, there are certain limitations when conducting tests using the traditional method of gravity-induced free fall. When the drop height is too large, the trajectory can be affected by air resistance, leading to deviations. On the other hand, when the drop height is too small, the rice grains may not bounce and instead roll on the surface, resulting in an inaccurate measurement of the restitution coefficient. To address these issues, this study proposed an impact collision method for determining the restitution coefficient of rice grains. It investigated the factors influencing the restitution coefficient, providing valuable insights for the accurate measurement of the restitution coefficient in rice grains.

Although this study accurately determined the relevant physical and contact parameters of rice grains and successfully calibrated the discrete element model, it did not incorporate the hollow characteristic of rice stems based on their actual physical features. Since this study did not specifically focus on the shape characteristics of straw stems, it has no significant impact on the research findings. However, if the discrete element modeling were to include the actual hollow feature of rice stems, it would increase the workload and computational time of the simulation, which would be unfavorable for the simulation process. Additionally, this study neglected the modeling of stem leaves, and the grain modeling was based on the Hertz–Mindlin (no slip) model, which cannot address the issues of grain impact fragmentation and residual characteristics during threshing. In future research, it would be beneficial to establish an accurate whole rice plant model using three-dimensional laser scanning and refine the modeling of rice grains to incorporate their fragmentation characteristics.

4. Conclusions

• The basic physical parameters, static friction coefficient, rolling friction coefficient, elastic modulus, and restitution coefficient of the whole rice plant were determined.
• A discrete element model of the whole flexible upright rice plant was established based on the Hertz–Mindlin (no slip) and Hertz–Mindlin with bonding models. A total of 43,214 particles were used for filling, resulting in the formation of 53,473 bonding contacts.
• Taking the calibration of contact parameters between rice grains and steel as an example, the calibration process was analyzed in detail, and the calibration error of the natural repose angle was determined to be 0.68%. Similarly, taking the calibration of contact parameters between the main stems and steel as an example, the calibration process was analyzed in detail, and the calibration errors of the natural repose angle and deflection were determined to be 2.76% and 2.33%, respectively.
• This study provides a fundamental understanding of the mechanical response of rice plants under mechanical coupling, thereby establishing a solid foundation for further research in this field. Moreover, the insights gained from this study can serve as a valuable reference for the development of discrete element models for plant structures in other crop species.