Calibration and Veriﬁcation Test of Lily Bulb Simulation Parameters Based on Discrete Element Method

: In the simulation analysis of the lily harvesting process, the intrinsic parameters of the lily bulb and the contact parameters between the lily bulb and the lily mechanized harvesting equipment (Q235 steel) are deﬁcient. Thus, the three-axis size, density, moisture content, Poisson’s ratio, elastic modulus, and other parameters of lily bulbs are measured in this paper with lily bulbs as the research object. Moreover, the discrete element model of the lily bulb was established using 3D scanning technology. The contact parameters between the lily bulb and Q235 steel were calibrated through bench test and simulation parameter test. The relative error between the measured value of the lily bulb accumulation angle and the simulated value is taken as response value to calibrate three parameters (collision recovery coefﬁcient, static friction coefﬁcient, and dynamic friction coefﬁcient between lily bulbs). A regression model of the relative error of the stacking angle and three parameters is established, and the response surface is optimized. The results demonstrate that collision recovery coefﬁcient, static friction coefﬁcient, and dynamic friction coefﬁcient between lily bulb and Q235 steel are 0.301, 0.423, and 0.063, respectively; these coefﬁcients between lily bulbs are 0.455, 0.425, and 0.158, respectively. Additionally, a better combination of parameters is adopted to perform the simulation stacking test. The measured stacking angle is 32.31 ◦ , which is 0.34% in error with the stacking angle obtained by the physical stacking test. The test results suggest that the discrete element model and contact parameters of the lily bulb can be used in the discrete element simulation test. Furthermore, these research results could provide references for simulation tests, such as mechanized harvesting and post-harvest processing, of lily. Author Contributions: Conceptualization, Z.D. and M.W.; methodology, Z.D.; software, Z.F.; valida-tion, Y.Q., Z.F. and Z.D.; formal analysis, M.W.; investigation, Y.Q.; resources, M.W.; data curation, Z.D.; writing—original draft preparation, Z.D.; writing—review and editing, Z.D.; visualization, Z.D.; supervision, M.W.; project administration, Z.F.; funding acquisition,


Introduction
Lily, belonging to the family Liliaceae, is of edible and medicinal significance, exhibiting high nourishing effect and economic value [1]. However, lilies are generally planted to a depth of 15-25 cm, causing soil disturbance during machine harvesting and increasing the chance of damage on lily bulbs [2]. A numerical simulation study should be conducted on the mechanical harvesting process of lily bulbs to achieve high-efficiency and lowdamage harvesting of lily bulbs and thus to provide a reference for the harvest mechanical design of lily bulbs. Before a numerical simulation study of the lily bulb is conducted, the intrinsic parameters of the lily bulb need to be determined [3][4][5]. At the same time, the contact parameters between the lily bulb and the machine (Q235) need to be calibrated. The parameter calibration can be completed through a virtual simulation test and physical test.
At present, researchers around the world have conducted a significant amount of research on the model establishment and parameter calibration of agricultural bulk materials. Certain progress has been made in discrete element modeling and parameter calibration of related bulk materials, mainly including material grains, potatoes, animal manure, and soil 2 of 18 particles [6][7][8][9][10][11]. Coetzee et al. [12] employed the shear test and the confined compression test to complete the calibration of the friction coefficient and stiffness coefficient of the corn grains while performing physical test verification. Ghodki et al. [13] established a soybean discrete element model and obtained the contact parameters of the model by comparing the characteristics of soybean shape and angle of repose. Ucgul et al. [14,15] achieved the calibration of soil cohesive force and non-cohesive soil simulation parameters by combining Hertz-Mindlin and Hysteretic spring contact models, and solved the problems of plastic deformation under stress. Hao Jianjun et al. [16] built a discrete element model of Ma yam based on the accumulation of bimodal distribution. Shi Ronglin et al. [17] designed a flax seed model based on the triaxial size of flax seeds and measured the physical parameters of three different flax seeds. Peng Caiwang et al. [18] simulated the pile-up angle of the pig manure organic fertilizer particles treated with Heishui Fly through the steepest climbing test and the Box-Behnken test. Lu Fangyuan et al. [19] calibrated the main contact parameters of the discrete elements of rice buds under different moisture contents depending on the experiment and simulation measurement of the friction angle of the rice buds. Research has demonstrated that it is feasible to obtain intrinsic parameters of bulk materials using the discrete element method. However, there is little research on the calibration of the discrete element model parameters for lily bulbs and machinery.
In sum, Longya lily planted in the Longhui area of Shaoyang, Hunan with a moisture content of 64.3% is used as a sample in this paper. Combined with bench tests and simulation tests, the establishment of the discrete element model of lily bulb was fulfilled using 3D reverse scanning to calibrate the contact parameters between the lily bulb and the machine (Q235 steel). Furthermore, the contact parameters between the lily bulbs were determined using a cylinder lifting, quadratic regression orthogonal rotation combination test, and response surface optimization methods. The pile-up angle test combining the physical test and the virtual test was conducted to verify the accuracy of lily bulb simulation parameters. There are two main innovations in this article. On the one hand, we established a discrete element model of lily bulb by using 3D reverse simulation engineering technology. On the other hand, we have verified the accuracy of the experimental values through a combination of virtual simulation experiments and physical experiments. Therefore, this research results can provide an accurate reference for the physical properties of lily in different stages, such as machine harvesting and postpartum processing.

Triaxial Size Determination of Lily Tubers
According to the requirements of discrete element simulation, the basic physical parameters and contact parameters of the lily tubers need to be provided. The basic physical parameters include the triaxial size, mass density, elastic modulus, Poisson's ratio, and moisture content of lily tubers. Contact mechanic parameters are composed of static friction coefficient, rolling friction coefficient, and collision recovery coefficient. Therefore, the triaxial dimensions of lily tubers should be determined first. In this study, the Longya lily planted in Longhui County, Shaoyang, Hunan Province was used as the research object. The five-point sampling method was used to sample Longya lily in Longhui County, Shaoyang to accurately establish a three-dimensional model of lily tubers. , 300 pieces were randomly selected, and a 111-101 type vernier caliper (precision 0.01 mm) was employed to measure their three-axis dimensions (maximum length L, maximum width W, maximum thickness T), the result being the averages collected. The three-axis size measurement of Longya Lily is illustrated in Figure 1. maximum thickness T), the result being the averages collected. The three-axis size m urement of Longya Lily is illustrated in Figure 1.

Moisture Content and Density
The moisture content of lily tubers was measured using the XF-110MA touchquick moisture meter (range 110 g). This instrument can be employed to automatic calculate the moisture content of materials according to their change quality. The w content of the lily tuber was measured after peeling off the lily tuber scales since the o all quality of a single lily tuber exceeded the maximum range of the instrument. Du measurement, lily scales were put into the instrument and heated until moisture con became constant, and then its value was collected. The measurement was repeated 5 ti to obtain the average value. The average moisture content of the lily tuber was calcul to be 64.3%.
DM-300, the density of lily tubers was measured using a density measuring in ment (range 0.01-300 g, accuracy 0.01 g). Similarly, a suitable number of lily scales taken for density determination. This instrument can be adopted to automatically ca late the density of lily tubers by measuring the mass of lily scales and volume of drain Similarly, the experiment was repeated 5 times to acquire the average value. The den of the lily tuber was calculated to be 986 Kg/m 3 .

Determination of Mechanical Characteristic Parameters of Lily
The stress state during the separation of fruit and soil should be analyzed to fur establish its discrete element model, so as to ensure that the lily tubers are not dama

Moisture Content and Density
The moisture content of lily tubers was measured using the XF-110MA touch-type quick moisture meter (range 110 g). This instrument can be employed to automatically calculate the moisture content of materials according to their change quality. The water content of the lily tuber was measured after peeling off the lily tuber scales since the overall quality of a single lily tuber exceeded the maximum range of the instrument. During measurement, lily scales were put into the instrument and heated until moisture content became constant, and then its value was collected. The measurement was repeated 5 times to obtain the average value. The average moisture content of the lily tuber was calculated to be 64.3%.
DM-300, the density of lily tubers was measured using a density measuring instrument (range 0.01-300 g, accuracy 0.01 g). Similarly, a suitable number of lily scales was taken for density determination. This instrument can be adopted to automatically calculate the density of lily tubers by measuring the mass of lily scales and volume of drainage. Similarly, the experiment was repeated 5 times to acquire the average value. The density of the lily tuber was calculated to be 986 kg/m 3 .

Determination of Mechanical Characteristic Parameters of Lily
The stress state during the separation of fruit and soil should be analyzed to further establish its discrete element model, so as to ensure that the lily tubers are not damaged during the separation of fruit and soil during harvest. The HY-0580 microcomputer electronic universal testing machine was employed to measure mechanical characteristic parameters of lily bulbs through bending, compression, and shear tests.

Poisson's Ratio
Lily tubers are too large in size and have thin and brittle skins and it is difficult to measure them directly with the general cylindrical sample compression method. Therefore, 10 lilies were randomly selected from the above samples, and lily scales with uniform texture were taken. Original dimensions of the lily scales in the width and the thickness directions were recorded. Hengyi HY-0580 universal material tensile and compression testing machine were taken to apply pressure (loading speed 0.1 mm/s) along with the width and thickness directions of the lily scale until it cracked. The universal material testing machine and an electronic vernier caliper were used to record the deformation in the width direction and the deformation in thickness of the lily at the limit of axial load cracking, respectively [20]. Then, Poisson's ratio of lily was calculated using Formula (1), and the average results were taken. The measured Poisson's ratio of lily was 0.426.
where µ denotes the Poisson's ratio, indicates the deformation in the width direction of the lily scales, mm; represents the deformation in the thickness direction of the lily scales, mm; W 1 and W 2 designate the width of the lily scales before and after loading, respectively, mm; T 1 and T 2 refer to the thickness of the lily scales before and after loading, respectively, mm.

Elastic Modulus and Shear Modulus
Elastic modulus is a scale used to measure the ability of a material to resist elastic deformation. The elastic modulus must be measured first to calculate the shear modulus. In this test, 10 lilies were selected with scales taken off and uniform texture. Afterward, the A111-101 vernier caliper was used to measure their thickness before compression and then placed freely on the circle of the Hengyi HY-0580 universal material tensile and compression testing machine. On the platform, a circular indenter with a diameter of 5 mm was used. The speed before and during compression was 0.04 mm/s and 0.02 mm/s, respectively. The force (F)-deformation (∆L) data was read and repeated for the selected 10 lilies. In the above test, the elastic and shear modulus of lily was calculated by Formulas (2) and (3): where E denotes the modulus of elasticity, MPa; F refers to the axial load on the lily, N; S is the contact area, mm 2 , the diameter of the circular indenter is 5 mm, and the contact area with the lily scale is 11.465 mm 2 ; ∆L indicates the lily scale subjected to the deformation after compression, mm; L represents the thickness of the lily scale before compression, mm; G is the shear modulus, MPa; designates the Poisson's ratio of the lily scale. Through 10 measurements, the average value of lily's elastic modulus was obtained to be 5.25 MPa, and the average value of the shear modulus is calculated to be 1.85 MPa.

Establishment of Discrete Element Model of Lily Bulb
A lily bulb is irregular in shape, and conventional modeling methods cannot accurately restore its true characteristics. In this paper, lily bulbs of moderate size (as shown in Figure 2a) are selected as the research object to accurately establish a three-dimensional model of lily tubers and improve the authenticity of the simulation experiment. Furthermore, the three-dimensional scanning technology was applied, and an SP01-3D threedimensional scanner (measurement accuracy 0.02 mm, scanning Range 10~1000 mm) was employed to scan the outer contour of the lily. Besides, 3D coordinates of the outer surface of the lily were accurately obtained to generate point cloud data and then exported to Geomagic Studio software (Geomagic Studio V2020, Rock Hill, NM, USA) for merging and splicing to obtain the lily model. Finally, the lily bulb model was imported into GOM Inspect (GOM Inspect, Braunschweig, Germany). The sharp and noisy points were sharpened using the software to obtain the three-dimensional lily model (Figure 2b) [21][22][23]. Next, the final lily 3D model was imported into EDEM2020 software (EDEM2020, Tory, USA) and filled with 149 spherical particles with radii ranging from 2-10 mm until the lily 3D model was tightly filled and there was no space to fill (Figure 2c). rately restore its true characteristics. In this paper, lily bulbs of moderate size (as shown in Figure 2a) are selected as the research object to accurately establish a three-dimensional model of lily tubers and improve the authenticity of the simulation experiment. Furthermore, the three-dimensional scanning technology was applied, and an SP01-3D three-dimensional scanner (measurement accuracy 0.02 mm, scanning Range 10~1000 mm) was employed to scan the outer contour of the lily. Besides, 3D coordinates of the outer surface of the lily were accurately obtained to generate point cloud data and then exported to Geomagic Studio software (Geomagic Studio V2020, Rock Hill, NM, USA) for merging and splicing to obtain the lily model. Finally, the lily bulb model was imported into GOM Inspect (GOM Inspect, Braunschweig, Germany) . The sharp and noisy points were sharpened using the software to obtain the three-dimensional lily model (Figure 2b) [21][22][23]. Next, the final lily 3D model was imported into EDEM2020 software (EDEM2020 , Tory, United States) and filled with 149 spherical particles with radii ranging from 2-10 mm until the lily 3D model was tightly filled and there was no space to fill (Figure 2c).

Calibration of Contact Parameters between Lily Bulb and Q235 Steel
In this paper, both bench and simulation tests are compared to calibrate the contact parameters of lily bulbs and Q235 steel and thus ensure the reference-ability of the discrete element simulation test. During the harvesting process, lily bulbs and Q235 steel material parts were in contact with each other. Additionally, the contact parameters between lily bulbs and Q235 steel and between lily bulbs are calibrated. Discrete element parameter calibration mainly adopted the collision bounce, inclined surface rolling, inclined surface slip, and accumulation tests. During the process of calibrating the contact parameters of lily bulbs and materials, the moisture content of lily bulbs had an essential influence on it. Therefore, lily bulbs with a moisture content of 64% after harvest were selected for the experiment. The simulation test parameters in EDEM software are exhibited in Table 2.

Calibration of Contact Parameters between Lily Bulb and Q235 Steel
In this paper, both bench and simulation tests are compared to calibrate the contact parameters of lily bulbs and Q235 steel and thus ensure the reference-ability of the discrete element simulation test. During the harvesting process, lily bulbs and Q235 steel material parts were in contact with each other. Additionally, the contact parameters between lily bulbs and Q235 steel and between lily bulbs are calibrated. Discrete element parameter calibration mainly adopted the collision bounce, inclined surface rolling, inclined surface slip, and accumulation tests. During the process of calibrating the contact parameters of lily bulbs and materials, the moisture content of lily bulbs had an essential influence on it. Therefore, lily bulbs with a moisture content of 64% after harvest were selected for the experiment. The simulation test parameters in EDEM software are exhibited in Table 2.

Collision Recovery Coefficient
Collision recovery coefficient is a parameter of the ability of an object to recover after contact collision deformation. The collision recovery coefficient between lily bulbs and Q235 steel plate was calibrated by the collision bounce test [25]. With the purpose of ensuring that lily bulbs do not break after being dropped and the bounce height is easy to distinguish, a preliminary test was conducted to reveal that the best height of lily bulbs to fall was 300 mm. As illustrated in Figure 3a, the Q235 steel sheet was placed horizontally, and the lily bulb was freely dropped from the height h 1 of 300 mm. The lily bulb collides with the Q235 steel plate, and the highest bounce height of the lily bulb was recorded by high-speed video. This was repeated 5 times to obtain the average value (highest bounce height h 2 = 27.5 mm).
Collision recovery coefficient is a parameter of the ability of an object to recover after contact collision deformation. The collision recovery coefficient between lily bulbs and Q235 steel plate was calibrated by the collision bounce test [25]. With the purpose of ensuring that lily bulbs do not break after being dropped and the bounce height is easy to distinguish, a preliminary test was conducted to reveal that the best height of lily bulbs to fall was 300 mm. As illustrated in Figure 3a, the Q235 steel sheet was placed horizontally, and the lily bulb was freely dropped from the height h1 of 300 mm. The lily bulb collides with the Q235 steel plate, and the highest bounce height of the lily bulb was recorded by high-speed video. This was repeated 5 times to obtain the average value (highest bounce height h2 = 27.5 mm). The static friction factor (x2) and rolling friction factor (x3) between lily bulbs and Q235 steel, as well as the collision recovery coefficient (X1) and static friction factor (X2) and rolling friction factor (X3) between lily bulbs, have no effect on the bounce height. The values of x2, x3, X1, X2, and X3 were all set to 0 to avoid interference in the EDEM simulation test. After pre-simulation tests, the range of the collision recovery coefficient x1 of the lily bulb and Q235 steel was determined to be 0.20~0.40, and the step length was 0.04. There were 6 sets of simulation tests. Each set of tests were repeated 5 times to obtain the average value. The test design scheme and results are presented in Table 3, in which y1 denotes the simulation value of the highest bounce height of Q235 steel plate. The curve fitting was performed on the test data in Table 3 to obtain the relationship between the collision recovery coefficient of lily bulb-Q235 steel and lily bulb-clay soil The static friction factor (x 2 ) and rolling friction factor (x 3 ) between lily bulbs and Q235 steel, as well as the collision recovery coefficient (X 1 ) and static friction factor (X 2 ) and rolling friction factor (X 3 ) between lily bulbs, have no effect on the bounce height. The values of x 2 , x 3 , X 1 , X 2 , and X 3 were all set to 0 to avoid interference in the EDEM simulation test. After pre-simulation tests, the range of the collision recovery coefficient x 1 of the lily bulb and Q235 steel was determined to be 0.20~0.40, and the step length was 0.04. There were 6 sets of simulation tests. Each set of tests were repeated 5 times to obtain the average value. The test design scheme and results are presented in Table 3, in which y 1 denotes the simulation value of the highest bounce height of Q235 steel plate. The curve fitting was performed on the test data in Table 3 to obtain the relationship between the collision recovery coefficient of lily bulb-Q235 steel and lily bulb-clay soil and the maximum bounce height in the simulation test. The second-degree polynomial fitting curve was acquired, as presented in Figure 4. The fitted curve equation (Equation (4)) is: (4) and the maximum bounce height in the simulation test. The second-degree poly fitting curve was acquired, as presented in Figure 4. The fitted curve equation (E (4)) is: y1 = 159.4x1 2 + 16.94x1 + 7.811 The determination system of Equation (4) was R 2 = 0.981, which is close to 1, re the high equation fitting reliability. The actual measured value of 27.5 mm of the lil Q235 steel plate was substituted into the Equation (4) to obtain x1 = 0.301, and inp the EDEM software for verification, repeated 5 times. Afterward, the average val calculated. According to the test, the largest rebound height value was 26.24 mm relative error value of 4.96%. The above comparison demonstrated that the sim value was basically consistent with the bench test value. In the EDEM simulation t coefficient of recovery between the lily bulb and Q235 steel was determined to 0.301.

Coefficient of Static Friction
The static friction factor is the ratio of the maximum static friction force expe by the object to the normal pressure, and it can be examined by the inclined plane m [26] The static friction coefficient between the lily bulb and Q235 steel can be calibr the slope sliding method. The comparative test is exhibited in Figure 5. In the beg Q235 steel plate was placed horizontally, and a lily bulb was placed on the Q2 plate. Thus, the Q235 steel plate rotated slowly and uniformly around one end. W lily bulb started to slip, it stopped rotating. The frame test was repeated 5 times to the average value. The measured value of the tilt angle of the inclined plate w 23.03°. The determination system of Equation (4) was R 2 = 0.981, which is close to 1, reflecting the high equation fitting reliability. The actual measured value of 27.5 mm of the lily bulb-Q235 steel plate was substituted into the Equation (4) to obtain x 1 = 0.301, and input into the EDEM software for verification, repeated 5 times. Afterward, the average value was calculated. According to the test, the largest rebound height value was 26.24 mm, with a relative error value of 4.96%. The above comparison demonstrated that the simulation value was basically consistent with the bench test value. In the EDEM simulation test, the coefficient of recovery between the lily bulb and Q235 steel was determined to be x 1 = 0.301.

Coefficient of Static Friction
The static friction factor is the ratio of the maximum static friction force experienced by the object to the normal pressure, and it can be examined by the inclined plane method. [26] The static friction coefficient between the lily bulb and Q235 steel can be calibrated by the slope sliding method. The comparative test is exhibited in Figure 5. In the beginning, Q235 steel plate was placed horizontally, and a lily bulb was placed on the Q235 steel plate. Thus, the Q235 steel plate rotated slowly and uniformly around one end. When the lily bulb started to slip, it stopped rotating. The frame test was repeated 5 times to acquire the average value. The measured value of the tilt angle of the inclined plate was ϕ 1 = 23.03 • .
The rolling friction factor between lily bulbs and Q235 steel (x 3 ), as well as the collision recovery coefficient (X 1 ) and static friction factor (X 2 ) and rolling friction factor (X 3 ) of lily bulbs, had no effect on the tilt angle of the inclined plate. In the EDEM simulation test, the values of x 3 , X 1 , X 2 , and X 3 were all set to 0 to avoid interference. The calibrated parameters were used in this research. The collision recovery coefficient between lily bulb and Q235 steel was x 1 = 0.301. After the pre-simulation test, the static friction coefficient x 2 between lily bulbs and Q235 steel in the range of 0.30~0.55, with a step length of 0.05, was used to conduct 6 sets of simulation tests. Each set of simulation tests were repeated 5 times to obtain the average value. The test design scheme and results are illustrated in the table, where y 2 denotes the simulation value of the inclination angle of Q235 steel plate. [26] The static friction coefficient between the lily bulb and Q235 steel can be calibrated by the slope sliding method. The comparative test is exhibited in Figure 5. In the beginning, Q235 steel plate was placed horizontally, and a lily bulb was placed on the Q235 steel plate. Thus, the Q235 steel plate rotated slowly and uniformly around one end. When the lily bulb started to slip, it stopped rotating. The frame test was repeated 5 times to acquire the average value. The measured value of the tilt angle of the inclined plate was φ1 = 23.03°. The curve fitting was performed on the test data in Table 4 to acquire the relationship between the static friction factor and the inclination angle of the lily bulb and Q235 steel during the simulation test. The second-degree polynomial fitting curve was presented in Figure 6. The curve equation (Equation (5)) is:  The rolling friction factor between lily bulbs and Q235 steel (x3), as well as th sion recovery coefficient (X1) and static friction factor (X2) and rolling friction facto lily bulbs, had no effect on the tilt angle of the inclined plate. In the EDEM simulati the values of x3, X1, X2, and X3 were all set to 0 to avoid interference. The calibrated eters were used in this research. The collision recovery coefficient between lily b Q235 steel was x1 = 0.301. After the pre-simulation test, the static friction coeffi between lily bulbs and Q235 steel in the range of 0.30~0.55, with a step length of 0. used to conduct 6 sets of simulation tests. Each set of simulation tests were rep times to obtain the average value. The test design scheme and results are illustrate table, where y2 denotes the simulation value of the inclination angle of Q235 steel The curve fitting was performed on the test data in Table 4 to acquire the relat between the static friction factor and the inclination angle of the lily bulb and Q2 during the simulation test. The second-degree polynomial fitting curve was prese Figure 6. The curve equation (Equation (5)) is: y2 = −72.643x2 2 + 106x2 − 8.8057  The coefficient of determination R 2 of Equation (5) is 0.9887, which is clo demonstrating the highly reliable fitting of the equation. Substituting the measure of the inclination angle of the bench test (23.03°) into Equation (5), x2 = 0.423 was o and then input into EDEM for verification. The test was repeated 5 times to ob average value. The simulated value of the inclination angle measured by the sim test was 22.05°; the error value between the simulated and measured values was The above comparison suggested that simulated and bench test values were the s The above comparison suggested that simulated and bench test values were the same. In the EDEM simulation test, the static friction factor x 2 between the lily bulb and Q235 steel was determined to be 0.423.

Rolling Friction Factor
The rolling friction coefficient is a vital contact parameter between the lily bulb and Q235 steel. The inclined rolling test is mainly performed to calibrate the rolling friction coefficient between the lily bulb and Q235 steel. This test is illustrated in Figure 7 [27]. This test was conducted by placing the lily bulb on an inclined surface with an inclination angle of ϕ 2 = 40 • , at a fixed height S 1 = 300 mm. Hence, the lily bulb rolls down the inclined surface with an initial speed of 0. The horizontal rolling distance of the lily bulb was measured when it rolled down and was completely still on the horizontal surface. The bench test was repeated 5 times to acquire the average value. The measured value of the horizontal rolling distance was S 2 = 149.8 mm.

Rolling Friction Factor
The rolling friction coefficient is a vital contact parameter between the lily bulb and Q235 steel. The inclined rolling test is mainly performed to calibrate the rolling friction coefficient between the lily bulb and Q235 steel. This test is illustrated in Figure 7 [27]. This test was conducted by placing the lily bulb on an inclined surface with an inclination angle of φ2 = 40°, at a fixed height S1 = 300 mm. Hence, the lily bulb rolls down the inclined surface with an initial speed of 0. The horizontal rolling distance of the lily bulb was measured when it rolled down and was completely still on the horizontal surface. The bench test was repeated 5 times to acquire the average value. The measured value of the horizontal rolling distance was S2 = 149.8 mm. The collision recovery coefficient between lily bulbs (X1), static friction (X2), and rolling friction (X3) had no effect on the horizontal rolling distance. In the EDEM simulation test, the values of X1, X2, and X3 were all set to 0 to avoid interference. The calibrated parameters were used in this research. The coefficient of recovery from the collision between lily bulb and Q235 steel was x1 = 0.301 while the value of static friction was x2 = 0.423. After the pre-simulation test, the rolling friction coefficient between lily bulb and Q235 steel x3 was in the range of 0.04~0.09, with a step length of 0.01. Moreover, 6 sets of the simulation were conducted with each set of experiments repeated 5 times to take the average value. The experimental design plan and results are provided in Table 5, where y3 represents the simulation value of the horizontal rolling distance.  Table 5 to obtain the relationship between the rolling friction factor and the horizontal rolling distance of the lily bulb and Q235 steel in the simulation test. The second-order polynomial fitting curve was obtained, as illustrated in Figure 8. The curve equation (Equation (6)) y3 is: The collision recovery coefficient between lily bulbs (X 1 ), static friction (X 2 ), and rolling friction (X 3 ) had no effect on the horizontal rolling distance. In the EDEM simulation test, the values of X 1 , X 2 , and X 3 were all set to 0 to avoid interference. The calibrated parameters were used in this research. The coefficient of recovery from the collision between lily bulb and Q235 steel was x 1 = 0.301 while the value of static friction was x 2 = 0.423. After the pre-simulation test, the rolling friction coefficient between lily bulb and Q235 steel x 3 was in the range of 0.04~0.09, with a step length of 0.01. Moreover, 6 sets of the simulation were conducted with each set of experiments repeated 5 times to take the average value. The experimental design plan and results are provided in Table 5, where y 3 represents the simulation value of the horizontal rolling distance. The curve fitting was performed on the test data in Table 5 to obtain the relationship between the rolling friction factor and the horizontal rolling distance of the lily bulb and Q235 steel in the simulation test. The second-order polynomial fitting curve was obtained, as illustrated in Figure 8. The curve equation (Equation (6)  The coefficient of determination of Equation (6) was R 2 = 0.9948, which was 1, indicating the high fitting reliability of the equation. By substituting the actua ured value of the horizontal rolling distance of 149.8 mm into Equation (6), x3 = 0.0 obtained. This value was further entered into EDEM for verification. The test was re 5 times to obtain the average values. The simulated value of the horizontal rolling d measured by the simulation test was 142.3 mm, with a measured relative error v 5.3%, suggesting that the calibration simulation results are basically consistent w bench test. In the EDEM simulation test, the coefficient of rolling friction between bulb and Q235 steel was determined to be x3 = 0.063.

Calibration of Contact Parameters between Lily Bulbs
Due to the extremely irregular dimensions of lily bulbs, it is difficult to determ contact parameters between lily bulbs. Measurement results of the related param lily bulbs obtained by traditional methods had large errors. In the process of free accumulation of lily bulbs, there were collision and friction forces between each li The collision recovery coefficient, static friction, and rolling friction factor betw bulbs all affected the shape of the pile angle. Therefore, the stacking angle test m [28] can be employed to calibrate the parameters in combination with the simulat bench tests. Simultaneously, the quadratic regression orthogonal rotation combin and the response surface optimization method were employed to obtain the bes simulation contact parameters of the lily bulbs.

Lily Bulb Stacking Test
The stacked bench test is shown in Figure 9a. The simulation test was cre EDEM, as indicated in Figure 9b. The test adopted Q235 steel material was used t the funnel. The upper and lower ends of the funnel had diameters of Φ450 mm an mm, respectively. A baffle was installed on the lower end of the funnel. A 1000 × 10 mm 3 bottom plate was placed at a distance of h3 = 400 mm from the bottom The coefficient of determination of Equation (6) was R 2 = 0.9948, which was close to 1, indicating the high fitting reliability of the equation. By substituting the actual measured value of the horizontal rolling distance of 149.8 mm into Equation (6), x 3 = 0.063 was obtained. This value was further entered into EDEM for verification. The test was repeated 5 times to obtain the average values. The simulated value of the horizontal rolling distance measured by the simulation test was 142.3 mm, with a measured relative error value of 5.3%, suggesting that the calibration simulation results are basically consistent with the bench test. In the EDEM simulation test, the coefficient of rolling friction between the lily bulb and Q235 steel was determined to be x 3 = 0.063.

Calibration of Contact Parameters between Lily Bulbs
Due to the extremely irregular dimensions of lily bulbs, it is difficult to determine the contact parameters between lily bulbs. Measurement results of the related parameters of lily bulbs obtained by traditional methods had large errors. In the process of free fall and accumulation of lily bulbs, there were collision and friction forces between each lily bulb. The collision recovery coefficient, static friction, and rolling friction factor between lily bulbs all affected the shape of the pile angle. Therefore, the stacking angle test method [28] can be employed to calibrate the parameters in combination with the simulation and bench tests. Simultaneously, the quadratic regression orthogonal rotation combined test and the response surface optimization method were employed to obtain the best value simulation contact parameters of the lily bulbs.

Lily Bulb Stacking Test
The stacked bench test is shown in Figure 9a. The simulation test was created in EDEM, as indicated in Figure 9b. The test adopted Q235 steel material was used to make the funnel. The upper and lower ends of the funnel had diameters of Φ450 mm and Φ200 mm, respectively. A baffle was installed on the lower end of the funnel. A 1000 × 1000 × 10 mm 3 bottom plate was placed at a distance of h 3 = 400 mm from the bottom of the funnel. The baffle and bottom plate materials were all Q235 steel. The lower end of the funnel was blocked with a baffle and filled with lily bulbs. The baffle quickly removed; then the lily bulbs fell freely, piled up on the bottom plate. When all the lily bulbs were still and the slope was stable, a camera was used to take a vertical photograph of the two sides of the lily bulb pile to obtain a picture of the pile angle, as illustrated in Figure 9a. Additionally, MATLAB was used to process the pile-up angle graph obtained from the experiment to reduce the manual measurement error [29]. The corresponding processing method was described as follows. First, grayscale processing was performed on the original image for binarization; afterward, the contour curve was extracted by canny edge detection and perform linear fitting; finally, the slope obtained by linear fitting was converted into an angle, that is, the accumulation angle of the physical accumulation test of lily bulbs. The contour extraction process is illustrated in Figure 10. The above test was repeated 5 times, and the average value was calculated to obtain the measured value of the accumulation angle of the lily bulb physical accumulation test ϕ 3 = 32.2 • .
Appl. Sci. 2021, 11, x FOR PEER REVIEW method was described as follows. First, grayscale processing was performed on th inal image for binarization; afterward, the contour curve was extracted by cann detection and perform linear fitting; finally, the slope obtained by linear fitting w verted into an angle, that is, the accumulation angle of the physical accumulation lily bulbs. The contour extraction process is illustrated in Figure 10. The above te repeated 5 times, and the average value was calculated to obtain the measured v the accumulation angle of the lily bulb physical accumulation test φ3 = 32.2°. In the EDEM simulation test, the calibrated contact parameters of the lily bul Q235 steel were used; the collision recovery coefficient X1, the static friction coeffic and dynamic friction coefficient X3 of the lily bulbs were selected as test factors. Th tive error of the stacking angle between the simulation and the bench tests was ta the test index. The calculation formula (Formula (7)) for the relative error of the st angle Δ was Appl. Sci. 2021, 11, x FOR PEER REVIEW 11 of 18 method was described as follows. First, grayscale processing was performed on the original image for binarization; afterward, the contour curve was extracted by canny edge detection and perform linear fitting; finally, the slope obtained by linear fitting was converted into an angle, that is, the accumulation angle of the physical accumulation test of lily bulbs. The contour extraction process is illustrated in Figure 10. The above test was repeated 5 times, and the average value was calculated to obtain the measured value of the accumulation angle of the lily bulb physical accumulation test φ3 = 32.2°.
(a) (b) In the EDEM simulation test, the calibrated contact parameters of the lily bulbs and Q235 steel were used; the collision recovery coefficient X1, the static friction coefficient X2, and dynamic friction coefficient X3 of the lily bulbs were selected as test factors. The relative error of the stacking angle between the simulation and the bench tests was taken as the test index. The calculation formula (Formula (7)) for the relative error of the stacking angle Δ was In the EDEM simulation test, the calibrated contact parameters of the lily bulbs and Q235 steel were used; the collision recovery coefficient X 1 , the static friction coefficient X 2 , and dynamic friction coefficient X 3 of the lily bulbs were selected as test factors. The relative error of the stacking angle between the simulation and the bench tests was taken as the test index. The calculation formula (Formula (7)) for the relative error of the stacking angle ∆ was where ϕ 3 denotes the simulated value of the accumulation angle, ( • ); ϕ 3 represents the measured value of the accumulation angle, ( • ).

Steepest Climb Test
The steepest climbing test was performed to determine the 0 level and optimal value interval of the quadratic regression orthogonal rotation combination test factors [30]. The steepest climbing test plan and results are presented in Table 6. The lily bulb stacking map is exhibited in Figure 11. Figure 11a-g displays the lily bulb stacking maps of the steepest climbing simulation test for a~g groups, respectively. Figure 11h provides the lily bulb stacking map of the bench test. Table 6. Scheme and results of steepest ascent experiment.

Experimental Factors
Test Results ϕ′ denotes the simulated value of the accumulation angle, (°); 3 ϕ represents measured value of the accumulation angle, (°).

Steepest Climb Test
The steepest climbing test was performed to determine the 0 level and optimal va interval of the quadratic regression orthogonal rotation combination test factors [30]. T steepest climbing test plan and results are presented in Table 6. The lily bulb stacking m is exhibited in Figure 11. Figure 11a-g displays the lily bulb stacking maps of the steep climbing simulation test for a~g groups, respectively. Figure 11h provides the lily b stacking map of the bench test. It can be observed from Table 6 that the relative error of the accumulation angle f decreases and then increases. The relative error of the simulation test of group d was smallest, and the optimal value interval was near the test of group d. Therefore, the t factors of groups c, d, and e are selected as the factors of −1, 0, and 1 level of the quadra regression orthogonal rotation combination test.

Quadratic Regression Orthogonal Rotation Combined Test
A three-factor quadratic regression orthogonal rotation combination test was co ducted to find the best parameter combination of collision recovery, static friction, a dynamic friction coefficients of lily bulbs in the EDEM simulation test. The simulation t factor codes are shown in Table 7. The test factor codes X1, X2, and X3 in Table 7 are t code values of bulb collision recovery, static friction, and dynamic friction coefficien respectively. The simulation test design scheme and results are offered in Table 7. The t result was the relative error Y between the simulated accumulation angle and the me ured accumulation angle. According to the binary regression fitting of the data in Table 8, the regression mod of the relative error Y of the accumulation angle and the collision recovery, the static fr tion, and the dynamic friction coefficients (X1), (X2) and (X3) of the lily bulbs were esta lished, expressed as Equation (8)   It can be observed from Table 6 that the relative error of the accumulation angle first decreases and then increases. The relative error of the simulation test of group d was the smallest, and the optimal value interval was near the test of group d. Therefore, the test factors of groups c, d, and e are selected as the factors of −1, 0, and 1 level of the quadratic regression orthogonal rotation combination test.

Quadratic Regression Orthogonal Rotation Combined Test
A three-factor quadratic regression orthogonal rotation combination test was conducted to find the best parameter combination of collision recovery, static friction, and dynamic friction coefficients of lily bulbs in the EDEM simulation test. The simulation test factor codes are shown in Table 7. The test factor codes X 1 , X 2 , and X 3 in Table 7 are the code values of bulb collision recovery, static friction, and dynamic friction coefficients, respectively. The simulation test design scheme and results are offered in Table 7. The test result was the relative error Y between the simulated accumulation angle and the measured accumulation angle. According to the binary regression fitting of the data in Table 8, the regression model of the relative error Y of the accumulation angle and the collision recovery, the static friction, and the dynamic friction coefficients (X 1 ), (X 2 ) and (X 3 ) of the lily bulbs were established, expressed as Equation (8): The analysis of variance in Table 9 suggested that the p-value of the regression model was less than 0.0001, the p-value of the lack-of-fit term was 0.1283, the coefficient of determination R 2 was 0.9223, the regression model was extremely significant, and the lackof-fit term was not significant. This reveals that there are no other main factors affecting the index. The coefficient of determination is close to 1, reflecting that the regression equation fits well. As demonstrated from Table 9, X 3 , X 2 X 3 , and X 3 2 have extremely significant effects on the accumulation angle of lily bulbs, and X 2 has the most significant effects on the accumulation angle of lily bulbs; this is due to p-value < 0.0001 and high F-value. The order of significance from large to small is: X 3 2 > X 3 > X 2 X 3 > X 2 . There was a quadratic nonlinear relationship between the test factors X 2 and X 3 and the index, and the interaction effect on the index was extremely significant. The response surface of the stacking angle is exhibited in Figure 12.

Parameter Optimization
The regression equation was solved using an optimization module of D 8.0.6 software (Design-Expert v8.0.6, Jiangsu, China) , with a minimum rel stacking angle as the goal. The response surface was analyzed, and our regr was optimized. The objective and constraint equations (Equation (9) Using the optimization function of Design-Expert 8.0.6, the relative erro accumulation angle obtained from the physical accumulation test of lily bu to perform the optimization, and 49 sets of optimized solutions were obtaine experiments were performed on the optimized solution. The simulation tes compared with the physical accumulation test results to find a set of optimi that are most similar to the physical stacking test stacking angle size and sha sion recovery coefficient between lily bulbs was 0.455, the static friction was dynamic friction coefficient was 0.158.

Stack Angle Verification Test
The optimized parameters were used for simulation verification (X1 0.427, X3 = 0.158). The optimization solution of this group was subjected t repeated simulation experiments. The average value was obtained and use stacking angle of 32.31° under this parameter combination. The stacking an tained from the physical stacking test was 0.34%. The comparison between t test and the physical test was illustrated in Figure 13. The results demonstra was no significant difference between the accumulation angle simulation a tion angle physical test results under the optimized simulation parameters. M

Parameter Optimization
The regression equation was solved using an optimization module of Design-Expert 8.0.6 software (Design-Expert v8.0.6, Jiangsu, China) , with a minimum relative error of stacking angle as the goal. The response surface was analyzed, and our regression model was optimized. The objective and constraint equations (Equation (9)) were: Using the optimization function of Design-Expert 8.0.6, the relative error value of the accumulation angle obtained from the physical accumulation test of lily bulbs was used to perform the optimization, and 49 sets of optimized solutions were obtained. Simulation experiments were performed on the optimized solution. The simulation test results were compared with the physical accumulation test results to find a set of optimized solutions that are most similar to the physical stacking test stacking angle size and shape. The collision recovery coefficient between lily bulbs was 0.455, the static friction was 0.427, and the dynamic friction coefficient was 0.158.

Stack Angle Verification Test
The optimized parameters were used for simulation verification (X 1 = 0.455, X 2 = 0.427, X 3 = 0.158). The optimization solution of this group was subjected to five sets of repeated simulation experiments. The average value was obtained and used to obtain a stacking angle of 32.31 • under this parameter combination. The stacking angle error obtained from the physical stacking test was 0.34%. The comparison between the simulation test and the physical test was illustrated in Figure 13. The results demonstrated that there was no significant difference between the accumulation angle simulation and accumulation angle physical test results under the optimized simulation parameters. Moreover, the shapes and angles of the stacking angles of the two are highly similar, indicating that the simulation parameters of this group are set accurately.

Conclusion
(1) Using 3D reverse simulation engineering technology, the 3D model of the lily bulb was obtained through 3D scanning reverse modeling. In the EDEM software, the automatic filling method was employed to obtain the discrete element model of the lily bulb composed of 149 lowest to highest 2~10 mm unequal diameter particles. The experiment verified the rationality of the model. (2) The numerical comparison between simulation and bench test was conducted to calibrate the contact parameters of the lily bulb and Q235 steel. The collision recovery coefficient between the lily bulb and Q235 steel was calibrated to be 0.301 by the impact bounce test method. Through the inclined plane sliding test method, the static friction coefficient of the lily bulb and the Q235 steel was obtained to be 0.423. With the inclined plane rolling test method, the rolling friction factor of the lily bulb and the Q235 steel was obtained as 0.063. Through experiments, Zengxi Li et al. [31] obtained the collision recovery coefficient, static friction coefficient and rolling friction coefficient of lily bulb and Q235 steel as 0.285, 0.437, and 0.069, respectively. The errors of the calibration parameters are 5.32%, 3.31%, and 9.52%. The comparison results show that the test results are true and reliable. (3) Through the accumulation test method, Design-Expert 8.0.6 software was used to process the test data. The response surface optimization method based on the quadratic regression orthogonal rotation combination test was employed to determine the best contact parameters of the lily bulbs in the EDEM simulation test. The collision recovery coefficient, static friction factor, and rolling friction factor of lily bulbs were 0.455, 0.425, and 0.158, respectively. In the accumulation angle simulation test, the accumulation angle value of the above-mentioned combination parameters was obtained as 32.31°; the relative error value between the accumulation angle simulation value and the physical accumulation angle was calculated to be 0.34%; the shape and angle of the accumulation angles of the two are highly similar. Through experiments, Zengxi Li et al. obtained the collision recovery coefficient, static friction coefficient and rolling friction coefficient of lily bulb as 0.438, 0.401, and 0.172, respectively. The errors of the calibration parameters are 3.74%, 5.65%, and 8.86%. The above verification test results demonstrated that the calibration results are true and reliable. This can provide a reference for the simulation of mechanized operations during the sowing and harvesting stage of lilies. (4) This article mainly provides simulation test parameters for the simulation test of the key equipment of the lily harvester. This research results can provide an accurate reference for the physical properties of lily in different stages, such as mechanical harvesting and post processing. In addition, it introduces the contact parameters between lily bulbs and Q235 steel. Similarly, according to the test method of the article, it can be applied to the test of contact parameters between other materials and metallic materials and non-metallic materials.

Conclusions
(1) Using 3D reverse simulation engineering technology, the 3D model of the lily bulb was obtained through 3D scanning reverse modeling. In the EDEM software, the automatic filling method was employed to obtain the discrete element model of the lily bulb composed of 149 lowest to highest 2~10 mm unequal diameter particles. The experiment verified the rationality of the model. (2) The numerical comparison between simulation and bench test was conducted to calibrate the contact parameters of the lily bulb and Q235 steel. The collision recovery coefficient between the lily bulb and Q235 steel was calibrated to be 0.301 by the impact bounce test method. Through the inclined plane sliding test method, the static friction coefficient of the lily bulb and the Q235 steel was obtained to be 0.423. With the inclined plane rolling test method, the rolling friction factor of the lily bulb and the Q235 steel was obtained as 0.063. Through experiments, Zengxi Li et al. [31] obtained the collision recovery coefficient, static friction coefficient and rolling friction coefficient of lily bulb and Q235 steel as 0.285, 0.437, and 0.069, respectively. The errors of the calibration parameters are 5.32%, 3.31%, and 9.52%. The comparison results show that the test results are true and reliable. (3) Through the accumulation test method, Design-Expert 8.0.6 software was used to process the test data. The response surface optimization method based on the quadratic regression orthogonal rotation combination test was employed to determine the best contact parameters of the lily bulbs in the EDEM simulation test. The collision recovery coefficient, static friction factor, and rolling friction factor of lily bulbs were 0.455, 0.425, and 0.158, respectively. In the accumulation angle simulation test, the accumulation angle value of the above-mentioned combination parameters was obtained as 32.31 • ; the relative error value between the accumulation angle simulation value and the physical accumulation angle was calculated to be 0.34%; the shape and angle of the accumulation angles of the two are highly similar. Through experiments, Zengxi Li et al. obtained the collision recovery coefficient, static friction coefficient and rolling friction coefficient of lily bulb as 0.438, 0.401, and 0.172, respectively. The errors of the calibration parameters are 3.74%, 5.65%, and 8.86%. The above verification test results demonstrated that the calibration results are true and reliable. This can provide a reference for the simulation of mechanized operations during the sowing and harvesting stage of lilies. (4) This article mainly provides simulation test parameters for the simulation test of the key equipment of the lily harvester. This research results can provide an accurate reference for the physical properties of lily in different stages, such as mechanical harvesting and post processing. In addition, it introduces the contact parameters between lily bulbs and Q235 steel. Similarly, according to the test method of the article, it can be applied to the test of contact parameters between other materials and metallic materials and non-metallic materials.